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Abstract 

A  sensitivity  analysis  of  empirical  parameters  used  in  physics-based  models  was 
completed  in  this  study  to  determine  their  effect  on  electron  densities  and  total 
electron  content  (TEC)  in  the  ionosphere.  The  model  used  was  the  Ionosphere- 
Plasmasphere  Model  (IPM)  developed  by  Utah  State  University.  The  empirical  pa¬ 
rameters  studied  include  the  0+/0  collision  frequency,  zonal  wind,  secondary  electron 
production,  nighttime  E  x  B  drifts,  and  tidal  structure.  The  sensitivity  analysis  was 
completed  by  comparing  a  default  run  of  the  IPM  to  a  run  with  the  parameter  ad¬ 
justed  for  three  geophysical  cases.  Many  of  the  comparisons  resulted  in  nonlinear 
changes  to  the  model  output.  Doubling  the  0+/0  collision  frequency  increased  NmF2 
up  to  250%  in  the  equatorial  anomalies.  Setting  the  zonal  winds  to  zero  resulted  in  a 
400%  increase  in  TEC  units  (TECU)  over  the  Southeast  Pacific  and  a  50%  decrease 
near  Madagascar.  It  was  found  that  changes  in  electron  density  and  TEC  are  directly 
proportional  to  how  daytime  production  is  scaled  to  account  for  secondary  electron 
production.  Decreasing  the  nighttime  downward  Ex  B  drift  resulted  in  a  160%-630% 
increase  in  low  latitude  TEC.  Finally,  modulating  the  E  x  B  drift  to  simulate  tidal 
forcing  reproduced  the  four- wave  pattern  of  enhanced  low  latitude  TEC. 
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SENSITIVITY  ANALYSIS  OF  EMPIRICAL  PARAMETERS  IN  THE 


IONOSPHERE-PLASMASPHERE  MODEL 

I.  Introduction 

1.1  Motivation  and  Background 

As  technology  increases,  the  number  of  satellites  and  communication  systems  that 
dwell  in  the  near-earth  environment  continue  to  increase  as  well.  Additionally,  com¬ 
munication  signals  both  to  and  from  space-based  platforms  also  propagate  through 
the  near-earth  environment.  These  satellites  and  communication  systems  are  directly 
affected  by  this  environment  known  as  the  ionosphere.  The  ionosphere  is  highly  vari¬ 
able  and  has  many  effects  on  these  systems  including  high  frequency  (HF)  radio  wave 
communication  blackouts,  Global  Positioning  System  (GPS)  degradation,  and  navi¬ 
gation  communication  impacts.  Today’s  increased  reliance  on  these  systems  requires 
more  accurate  ionospheric  forecasts  in  order  to  mitigate  these  effects. 

The  ionosphere  is  the  region  in  the  earth’s  atmosphere  that  consists  of  a  quasi¬ 
neutral  ionized  plasma.  The  structure  and  characteristics  of  the  ionosphere  are  highly 
dependent  on  solar  cycle,  season,  time  of  day,  geographic  location,  and  geomagnetic 
activity.  Models  are  used  to  forecast  how  these  variables  change  the  ionosphere. 
Ionospheric  model  types  include  empirical,  parameterized,  numerical,  tomographic, 
and  physics-based  models  ( Schunk  et  al,  2002). 

Physics-based  models  typically  solve  the  continuity,  momentum,  and  energy  equa¬ 
tions  for  ions  and  electrons  taking  into  account  all  of  the  chemical  and  transport  pro¬ 
cesses  thought  to  be  important.  If  the  equations  describing  the  physics  and  chemistry 
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of  the  ionosphere  are  correct,  the  model  output  should  describe  the  real  ionosphere. 
However,  the  accuracy  of  the  model  output  depends  on  the  quality  of  the  input. 
The  physics-based  models  also  require  atmospheric  and  magnetospheric  inputs  as 
boundary  conditions  and  feedback  mechanisms  due  to  the  inherent  coupling  of  these 
systems.  When  sufficient  data  is  not  available,  empirical  models  are  often  used. 

Empirical  models  use  a  statistical  solution  based  on  limited  physics  and  many 
years  of  observed  data  (climatology).  In  general,  the  empirical  models  are  more 
correct  than  the  physics-based  models  because  the  output  is  based  on  averaging  real 
data  instead  of  solving  simplified  physical  equations  ( Schunk ,  2010).  However,  when 
these  empirical  models  are  combined  with  the  physics-based  models,  erroneous  and 
inconsistent  features  can  occur  in  the  model  output  due  to  the  combination  of  the 
averaged  data  from  the  empirical  model  and  the  more  robust  physics  from  the  physics- 
based  model.  These  inconsistencies  need  to  be  corrected  via  data  assimilation  within 
the  model  or  by  adjusting  certain  physical  parameters  in  the  model  as  determined  by 
comparison  with  observed  data  prior  to  use  of  the  model. 

For  example,  when  the  physics-based  Ionosphere  Forecast  Model  (IFM)  was  com¬ 
pared  to  actual  data  as  part  of  the  validation  process,  it  was  found  that  certain 
model  parameters  had  to  be  adjusted  for  the  model  output  to  match  observations 
(Appendix  A).  The  parameters  that  needed  to  be  adjusted  included  the  0+/0  colli¬ 
sion  frequency,  nighttime  drifts,  daytime  production,  zonal  winds,  tidal  forcing,  and 
electron  densities  at  high  latitudes  ( Schunk ,  2010).  The  IFM  is  part  of  the  Global 
Assimilation  of  Ionospheric  Measurements  (GAIM)  project  developed  at  Utah  State 
University  (USU)  and  is  used  by  the  Air  Force  as  the  operational  model  to  specify 
the  ionosphere.  The  IFM  has  been  validated  against  observations  for  different  pa¬ 
rameters  by  Scherliess  et  al.  (2006),  Zhu  et  al.  (2006),  Decker  and  McNamara  (2007), 
Thompson  et  al.  (2006),  and  others. 
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1.2  Research  Objective 


The  objective  of  this  research  is  to  complete  a  sensitivity  analysis  of  empirical 
parameters  in  a  physics-based  model.  This  is  done  by  examining  how  adjusting 
certain  physical  parameters  affects  the  model  output.  The  model  used  in  this  project 
is  the  Ionosphere-Plasmasphere  Model  (IPM).  The  IPM  is  a  physics-based  ionospheric 
forecast  model  under  development  at  USU  for  use  in  the  GAIM  Full  Physics  model, 
which  is  slated  to  be  the  next  operational  ionospheric  forecast  model  for  the  Air 
Force.  The  physical  parameters  that  are  examined  are  the  0+/0  collision  frequency, 
nighttime  drifts,  daytime  production,  zonal  winds,  and  tidal  forcing.  The  default 
model  output  is  compared  to  adjusted  output  to  determine  the  extent  of  the  changes 
in  model  output  due  to  the  adjustment  of  each  physical  parameter.  This  comparison 
will  give  an  understanding  of  the  importance  of  correctly  specifying  these  physical 
parameters  in  the  physics-based  models. 

1.3  Document  Structure 

This  document  is  organized  into  five  chapters.  Chapter  II  provides  background 
information  on  the  creation  and  transport  of  the  ionosphere  and  a  detailed  discussion 
of  the  IPM.  Chapter  III  details  the  procedures  used  to  conduct  this  research  and 
outlines  the  selection  of  geophysical  conditions,  running  the  model,  and  data  compar¬ 
ison  techniques.  Chapter  IV  gives  the  results  and  analysis  of  adjusting  the  physical 
parameters  in  a  physics-based  model  while  Chapter  V  summarizes  these  results,  gives 
final  conclusions,  and  outlines  recommendations  for  future  work. 
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II.  Background 


2.1  Chapter  Overview 

This  chapter  provides  background  information  necessary  to  understand  this  project 
and  consists  of  two  sections.  The  first  section  describes  the  composition  of  the  iono¬ 
sphere  and  plasmasphere.  Both  formation  and  transport  of  the  ionosphere  will  be 
discussed  in  this  section.  The  second  section  presents  a  detailed  explanation  of  the 
IPM  including  the  domain,  resolution,  and  model  physics. 

2.2  The  Ionosphere  and  Plasmasphere 

The  ionosphere  is  the  region  in  the  earth’s  atmosphere  that  consists  of  a  quasi¬ 
neutral  ionized  plasma.  It  extends  from  about  50 km  to  beyond  2000/cm  ( Schunk  and 
Nagy,  2009).  In  the  ionosphere,  the  plasma  density  is  high  enough  to  reflect  radio 
frequencies  while  the  neutral  density  remains  high  enough  to  influence  plasma  motion 
and  chemistry.  The  ionosphere  is  broken  up  into  different  regions  based  on  local  peaks 
in  the  electron  density.  These  regions  include  the  D  region  from  50-90/cm,  the  E  region 
from  90-140 km,  the  Fi  region  from  140-200/cm,  and  the  F2  region  from  200-300/cm. 
The  region  above  the  F2  electron  density  peak  is  known  as  the  topside  ionosphere 
and  extends  from  300-1500/cm  at  mid-latitudes  (Figure  1).  The  plasmasphere  is 
the  region  between  the  ionosphere  and  the  magnetosphere  consisting  of  low  energy, 
relatively  dense  plasma  that  begins  where  H+  becomes  the  dominant  ion  and  ends 
at  3-6  earth  radii,  depending  on  geomagnetic  activity  ( Tascione ,  1994).  In  general, 
the  plasmasphere  lies  beneath  the  magnetic  field  line  that  maps  to  ±60°  magnetic 
latitude  at  the  surface.  The  structure  and  characteristics  of  the  ionospheric  regions 
are  highly  dependent  on  solar  cycle,  season,  time  of  day,  geographic  location,  and 
geomagnetic  activity. 
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Figure  1.  Vertical  structure  of  electron  density  and  dominant  ions  showing  the  iono¬ 
spheric  regions  and  the  plasmasphere  (Adapted  from  Schunk  and  Nagy  (2009)) 


The  production,  loss,  and  transport  of  ions  control  the  ion  and  electron  density 
distributions  in  the  ionosphere.  These  three  processes  can  be  described  and  modeled 
using  the  continuity  equation  (Equation  1).  This  mathematical  expression  equates 
the  time  rate  of  change  of  the  density  (^r)  to  the  transport  (—V  •  Nsus )  and  the 
production  (Ps)  and  loss  ( Ls )  of  a  species  s 


dNs 

dt 


— V  •  Nsus  +  (Ps  —  Ls ) 


(1) 


where  Ns  is  the  density  of  the  species  and  us  is  the  velocity  of  the  species.  In  general, 
transport  is  negligible  below  200 km  and  the  change  in  ion  concentration  is  described 
only  by  production  and  loss  mechanisms.  This  photochemical  equilibrium  is  valid  in 
the  D,  E,  and  Fi  regions.  Above  200 km  (F2  region),  production  and  loss  no  longer 
dominate  the  continuity  equation  and  transport  must  be  included.  In  the  topside 
ionosphere  and  plasmasphere,  production  and  loss  are  negligible  and  transport  is  the 
dominant  process  that  controls  ion  and  electron  concentrations. 
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Where  photochemical  equilibrium  applies  in  the  D,  E,  and  Fi  regions,  the  iono¬ 


sphere  can  be  represented  by  a  simplified  continuity  equation  as 


dNs 

dt 


Ps  —  Ls 


(2) 


where  the  change  in  species  concentration  is  dependent  only  on  production  and  loss 
processes.  Photoionization  is  the  primary  production  mechanism  in  these  regions 
with  chemical  reactions,  secondary  electron  production,  and  impact  ionization  acting 
as  additional  production  mechanisms. 

Photoionization  creates  ion-electron  pairs  by  the  absorption  of  solar  extreme  ul¬ 
traviolet  (EUV)  and  X-ray  radiation  by  the  neutrals.  EUV  and  X-ray  radiation 
encompasses  a  wavelength  range  from  approximately  bum  to  llOnm  and  are  part 
of  the  solar  flux  emitted  from  the  sun  (Figure  2).  While  EUV  and  X-ray  radiation 
account  for  only  a  small  portion  of  the  overall  solar  flux,  it  is  these  wavelengths  that 
vary  the  most  with  solar  cycle — over  a  100%  increase  in  irradiance  from  solar  mini- 


Figure  2.  Solar  irradiance  and  variability.  Solar  spectrum  from  lOnra  to  105nra  (black) 
as  compared  to  a  5770  K  blackbody  (red)  and  the  solar  spectrum’s  variability  between 
solar  maximum  and  solar  minimum  (blue,  right  axis).  EUV  radiation  is  much  more 
variable  with  solar  cycle  than  visible  or  infrared  radiation  (Adapted  from  Lean  (1991)) 
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mum  to  solar  maximum  (Lean,  1991) — making  the  structure  of  the  ionosphere  highly 
dependent  on  the  solar  cycle  as  well. 

Production  due  to  photoionization  also  depends  on  the  optical  depth  of  the  at¬ 
mosphere,  neutral  concentrations,  and  ionization  cross  sections.  The  photoionization 
production  rate  for  a  particular  species  s  is  calculated  using 

r^si 

Ps(z,x)  =  Ns(z)  I00(X)exp[-T(z,x,X)]crl(X)dX  (3) 

do 

where 


/00(A) 

is  the  solar  flux  at  the  top  of  the  atmosphere 

Ns(z) 

is  the  number  density 

(z,X,  A) 

is  the  optical  depth 

as(A) 

is  the  ionization  cross  section 

X 

is  the  solar  zenith  angle 

A  si 

is  the  threshold  wavelength  for  ionization 

The  optical  depth  determines  how  much  of  the  incoming  solar  radiation  passes 
through  the  atmosphere.  It  is  dependent  on  the  neutral  concentrations  as  well  as  the 
relevant  absorption  cross  sections  and  the  angle  of  incoming  solar  radiation  (the  solar 
zenith  angle).  The  highest  concentration  of  neutrals  occurs  in  the  lower  atmosphere 
with  the  density  of  each  species  falling  off  exponentially  as  a  function  of  mass,  i.e.  the 
heavier  species  (Ar,  O2,  N2,  and  O)  fall  off  faster  than  the  lighter  species  (H  and  He) 
(see  Appendix  B).  Because  of  this  neutral  structure,  the  optical  depth  is  largest  at 
low  altitudes  and  the  intensity  of  the  incoming  photon  flux  decreases  as  the  radiation 
penetrates  into  the  atmosphere.  The  combination  of  decreasing  neutral  densities  and 
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increasing  solar  flux  with  altitude  creates  a  region  of  maximum  ionization  known  as 
a  Chapman  layer  (Figure  3). 


Figure  3.  Chapman  layer  and  ionization  curve  (dashed  line).  Decreasing  neutral  den¬ 
sities  (<t i  ■  n(z))  and  increasing  solar  irradiance  (I(z))  with  altitude  create  an  ionization 
peak  known  as  a  Chapman  layer  at  F  region  altitudes 

This  region  of  maximum  ionization  is  also  highly  dependent  on  the  solar  zenith 
angle  (dependence  is  carried  in  the  optical  depth  term).  As  the  zenith  angle  increases, 
the  peak  ionization  rate  decreases  and  the  height  of  the  peak  increases.  (S chunk  and 
Nagy ,  2009).  Therefore,  local  noon  will  have  the  highest  ionization  rate  with  the 
ionization  peak  at  its  lowest  altitude;  as  the  sun  nears  the  horizon,  the  ionization  rate 
decreases  and  the  height  of  the  peak  increases.  Photoionization  is  reduced  significantly 
at  night  due  to  the  loss  of  solar  radiation.  Ionization  cross-sections  are  the  probability 
that  ionization  will  occur  when  an  atom  or  molecule  interacts  with  a  photon.  These 
cross-sections  are  dependent  on  the  photon  wavelength  and  the  neutral  species  that 
the  photon  is  interacting  with. 

Two  important  chemical  reactions  that  create  ions  include  the  charge  exchange 
reaction  between  0+  and  H  and  the  bimolecular  ion-molecule  reaction  between  0+ 


and  N2.  These  two  reactions  are  given  in  Equations  4  and  5  below. 


0+  +  H  < 

>  H+  +  0 

(4) 

o+  +  n2  - 

NO+  +  N 

(5) 

Secondary  electron  production  occurs  when  an  EUV  or  X-ray  photon  photoionizes 
an  atom  or  molecule  and  the  resulting  electron  has  enough  energy  to  then  ionize  a 
second  atom  or  molecule.  This  process  is  shown  in  Equation  6 

O  +  hv  — »  0+  +  e  (6) 

where  e  denotes  the  resulting  high  energy  electron.  The  higher  the  energy  of  the 
electron,  the  more  times  it  is  able  to  ionize  another  atom  or  molecule  as  it  moves 
down  through  the  ionosphere.  Figure  4  shows  production  rates  both  with  and  without 
secondary  electrons.  As  seen  in  this  figure,  secondary  electron  production  increases 
production  rates  throughout  the  ionosphere,  but  does  not  significantly  alter  the  height 
of  the  peak  density.  The  final  production  mechanism,  impact  ionization,  is  a  result 
of  particle  precipitation  and  occurs  primarily  at  high  latitudes. 

The  primary  loss  processes  in  the  D,  E,  and  Fi  regions  are  chemical  reactions  in¬ 
cluding  ion-ion  recombination,  radiative  recombination,  dissociative  recombination, 
ion- atom  interchange,  and  associative  detachment.  The  dominant  loss  mechanism  in 
a  region  is  determined  by  the  chemical  reaction  rates;  dissociative  recombination  dom¬ 
inates  in  the  E  region  and  ion-atom  interchange  dominates  in  the  Fi  region  ( S chunk 
and  Nagy,  2009).  Appendix  B  gives  more  details  on  reaction  rates  as  well  as  produc¬ 
tion  and  loss  mechanisms  in  the  ionosphere. 

Chemical  reactions  can  be  a  production  mechanism  for  one  species  while  at  the 
same  time  a  loss  mechanism  for  another  species.  This  results  in  the  continuity  equa- 


9 


400 


Figure  4.  Production  rates  with  ( - )  and  without  ( - )  secondary  electrons  for  0+ 

(black),  ^2+  (red),  and  ()2  (blue)  during  solar  minimum.  Secondary  electrons  increase 
the  production  rate  at  all  altitudes,  but  do  not  affect  the  altitude  of  the  peak  (Adapted 
from  Schunk  and  Nagy  (2009)) 


tion  for  one  species  being  coupled  to  many  other  continuity  equations.  Therefore,  a 
numerical  solution  is  needed  to  solve  even  these  simplified  continuity  equations. 

In  the  F2  region,  transport  of  the  plasma  becomes  important.  This  is  modeled 
using  the  full  continuity  equation  (Equation  1)  which  includes  the  transport  term 
(—’V-Nsus)  and  the  transport  velocity  us.  To  obtain  the  transport  velocity,  a  solution 
to  each  species’  momentum  equation  is  needed.  When  thermal  diffusion,  Coriolis 
force,  and  centripetal  force  are  all  assumed  to  be  negligible,  the  general  momentum 
equation  is  given  as 


Ps 


Dsus 

~dT 


=  — Vps  -  V  •  ts  +  psg  +  nsqs 


E  +  us  x  B 


+!>-(*-*-)  (7) 


where  the  subscripts  denote  the  species  s  and  the  target  t,  ps  is  the  mass  density, 
Ds/Dt  is  the  convective  derivative,  ps  is  the  pressure,  kb  is  Boltzmann’s  constant,  ts 
is  the  stress,  and  ust  is  the  collision  frequency. 
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To  understand  plasma  transport,  the  motion  can  be  simplified  into  two  cases; 
motion  parallel  to  (along)  the  magnetic  field,  and  motion  perpendicular  to  (across) 
the  magnetic  field.  The  overall  plasma  motion  is  a  combination  of  these  two  motions 
and  can  be  understood  by  looking  at  the  approximations,  simplifying  assumptions, 
and  solutions  to  Equation  7. 

For  motion  along  the  magnetic  field,  the  diffusion  approximation  (steady  state  and 
subsonic),  charge  neutrality  (ne  =  ni),  and  zero  current  ( neue  =  nffii)  assumptions 
are  made,  resulting  in  the  ambipolar  diffusion  equation 


*11 


-Da 


V||Uj  +  V||  Tp 

'H  -Lp 


mk  (v  ~  ffh 

2  kbTp  2nikbTp 


rrtj 

2khTp 


VjnXLr 


(8) 


where  the  ambipolar  diffusion  coefficient  (. Da )  and  plasma  temperature  (Tp)  are 

2  kbTp 

Te  +  Tj 
2 


Da  = 
TP  = 


Equation  8  shows  that  the  ions  move  along  the  magnetic  field  with  the  neutral  wind 
subject  to  density  gradients,  temperature  gradients,  gravity,  and  stress  forces. 

The  two  dominant  forces  in  this  equation  are  the  neutral  winds  and  the  density 
gradient.  The  parallel  density  gradient  (V||7i,t)  points  from  low  to  high  densities.  For 
0+  in  the  F2  region,  the  gradient  is  in  the  upward  direction  because  0+  is  not  present 
in  the  lower  portions  of  the  ionosphere.  However,  the  diffusion  coefficient  is  proceeded 
by  a  negative  and  the  resulting  motion  due  to  diffusion  is  down  the  field  lines.  Above 
the  F2  peak,  the  parallel  0+  density  gradient  is  downward  during  the  day  due  to  the 
ionization  peak  and  upward  at  night  due  to  the  decay  of  0+ .  The  resulting  flow  due 
to  diffusion  in  the  topside  ionosphere  is  then  upward  during  the  day  and  downward 
at  night  ( Schunk  and  Nagy ,  2009). 
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The  thermospheric  neutral  winds  at  F2  region  altitudes  are  primarily  controlled  by 
solar  heating.  In  general,  the  thermospheric  winds  are  directed  radially  away  from  the 
subsolar  point  with  daytime  meridional  winds  blowing  toward  the  poles  and  nighttime 
meridional  winds  blowing  toward  the  equator.  As  seen  in  Figure  5,  the  winds  are  weak 
near  the  subsolar  point  and  increase  in  magnitude  as  the  distance  from  the  subsolar 
point  increases.  Weak  daytime  winds  are  due  to  higher  ion  concentrations  resulting 
in  increased  ion  drag  on  the  neutral  wind.  Because  the  plasma  is  constrained  to 
follow  the  held  lines,  the  neutral  winds  push  the  plasma  down  the  held  lines  to  lower 
altitudes  during  the  day  and  up  the  held  lines  to  higher  altitudes  at  night. 


Figure  5.  Thermospheric  neutral  winds  at  300/, 7//  calculated  using  HWM  93.  Winds 
flow  radially  away  from  the  subsolar  point  and  are  weaker  on  the  day  side  (Adapted 
from  Prolss  (2004)) 


Motion  across  the  magnetic  held  can  be  looked  at  in  a  manner  similar  to  motion 
along  the  magnetic  held.  After  transforming  to  a  reference  frame  moving  with  the 
neutral  wind  and  separating  the  perpendicular  motion  into  its  two  possible  orienta¬ 
tions,  the  motion  across  the  magnetic  held  is  described  by 


1 


Us±  = 


1  +  K2 


—Ds_  1  _  W 

Vj_ps  H  g±  +  /rs-E_|_ 

Ps  V sn 


+ 


1 


1  +  l/K< 


[lip  +  Ug  +  Up]  (9) 
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where 


K  =  — —  is  the  gyrofrequency  to  collision  frequency  ratio 

V sn 

Ds  =  —  -  is  the  diffusion  coefficient 

/j,s  =  — — —  is  the  mobility 

and  where 

(10) 

(11) 

(12) 


-1 


Up  — 


Un 


Up  — 


(v_l  ps 


x  B 


nsqs 

B 2 

ms  ( 

g±  x  b) 

<ls 

B 2 

x  B) 

__  7  is 

B 2 


is  the  gradient  drift 
is  the  gravitational  drift 
is  the  electromagnetic  drift 


At  high  altitudes  where  collisions  are  negligible,  K  — >  oo  and  the  first  term  in 
Equation  9  drops  out  leaving  only  motion  perpendicular  to  both  the  magnetic  field  and 
the  force.  Electrons  and  ions  drift  across  the  magnetic  field  in  opposite  directions  in 
the  presence  of  pressure  gradients  and  gravity,  but  they  drift  together  in  the  presence 
of  a  perpendicular  electric  field.  Typically  the  drifts  due  to  pressure  gradients  and 
gravity  are  small  and  the  dominant  drift  is  the  electromagnetic  drift.  At  low  latitudes, 
the  electromagnetic  drift  is  caused  by  the  dynamo  electric  field  created  in  the  E  region 
of  the  ionosphere  and  transmitted  along  the  highly  conductive  dipole  field  lines  to 
the  F  region.  This  electric  field  is  generated  by  thermospheric  winds  as  the  ions 
are  dragged  across  the  magnetic  field.  The  dynamo  electric  field  is  eastward  during 
the  day  with  a  resulting  upward  E  x  B  drift;  the  opposite  occurs  at  night  with  a 
westward  electric  field  and  a  downward  drift  (S chunk  and  Nagy ,  2009).  The  daytime 
combination  of  upward  drifts  and  downward  diffusion  is  responsible  for  the  ionization 
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peaks  on  both  sides  of  the  geomagnetic  equator  known  as  the  equatorial  anomaly 
(Figure  6). 


niffncir»r»  - -  niffiicinn 


Figure  6.  Formation  of  the  equatorial  anomaly  (pink)  is  a  result  of  upward  E  x  B  drifts 
and  diffusion  down  the  magnetic  field  lines 

At  low  altitudes  where  the  magnetic  field  strength  is  small  compared  to  the  colli¬ 
sion  frequency  K  — y  0,  and  the  second  term  in  Equation  9  becomes  negligible.  This 
occurs  because  of  the  high  concentration  of  neutrals  in  the  lower  atmosphere.  The 
resulting  plasma  motion  is  perpendicular  to  the  magnetic  field  but  in  the  same  di¬ 
rection  as  the  force  causing  the  motion.  At  intermediate  locations  where  neither  of 
the  above  approximations  is  valid  (i.e.  K  ~  1),  the  plasma  motion  is  a  combination 
of  the  motion  perpendicular  to  the  magnetic  field  and  parallel  to  the  force  ( Kelley , 
1989).  Plasma  movement  occurs  throughout  the  entire  ionosphere,  but  becomes  sig¬ 
nificant  in  the  F2  region  and  is  the  dominant  process  in  the  topside  ionosphere  and 
plasmasphere.  The  full  derivation  of  the  momentum  equation  solutions  are  included 
in  Appendix  B. 

Transport  and  diffusion  are  important  in  the  F2  region  because  this  is  where 
the  transition  from  chemical  equilibrium  in  the  D,  E,  and  Fi  regions  to  diffusive 


14 


equilibrium  in  the  topside  ionosphere  occurs.  The  location  and  magnitude  of  the  F2 
peak  is  dependent  on  where  the  chemical  and  diffusive  equilibrium  solutions  are  equal 
(Figure  7).  Changes  to  either  solution  by  changing  such  factors  as  the  the  diffusion 
rate  or  the  neutral  densities  can  change  the  magnitude  and  altitude  of  the  F2  peak. 


Figure  7.  Chemical  versus  diffusive  equilibrium  in  the  F2  region.  The  0+  density  peak 
occurs  where  chemical  and  diffusive  equilibrium  are  equal  (Adapted  from  Schunk  and 
Nagy  (2009)) 


Where  transport  and  diffusion  are  the  dominant  processes  in  the  topside  iono¬ 
sphere  and  plasmasphere,  the  ion  and  electron  concentrations  can  be  described  by  a 
diffusive  equilibrium  approximation.  To  simplify  this  approximation,  Equation  8  is 
transformed  from  a  Cartesian  coordinate  system  to  coordinates  along  the  magnetic 
field  as 

1  drij  _  rrtjg  1  dTp  dr^/dr  (un  -  fq)  .  . 

rii  dr  2 kfTp  Tp  dr  2 riikfTp  Da 

Because  the  diffusion  coefficient  is  inversely  proportional  to  the  neutral  density,  and 
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the  neutral  density  decreases  exponentially  with  altitude,  the  last  term  in  Equation  13 
becomes  negligible  as  altitude  increases  (Da  — »  oo).  At  high  altitudes  the  stress 
term  is  small  and  can  also  be  neglected  (t;||  — »  0).  With  these  approximations,  the 
momentum  equation  becomes  the  classical  diffusive  equilibrium  equation 

1  dni 
^  dr 

where  Hp  is  the  plasma  scale  height 


For  an  isothermal  atmosphere  and  ignoring  the  variation  of  gravity  with  altitude, 
Equation  14  is  integrated  to  give 


1  1  dTp 

1TP  ~  Tp~dr 


2  kbTp 


(14) 


Uj(r)  =  nj(r0)  e  ^  r°l/Hp 


(15) 


where  r  indicates  an  altitude  dependence  and  the  subscript  0  corresponds  to  a  ref¬ 
erence  altitude,  typically  the  height  of  the  F2  peak.  Equation  15  shows  that  in  this 
approximation  the  major  ion  density  decreases  exponentially  with  altitude  at  a  rate 
governed  by  the  plasma  scale  height  ( Schunk  and  Nagy ,  2009). 

In  addition  to  the  continuity  and  momentum  equations,  the  temperature  of  each 
species  due  to  heat  transport  is  important  in  describing  the  true  nature  of  the  iono¬ 
sphere.  This  temperature  change  is  described  by  the  energy  equation 
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(16) 


where  qs  is  the  heat  flow,  Qs  and  Ls  are  the  local  heating  and  cooling  rates,  respec¬ 
tively,  and  all  other  variable  are  as  described  previously.  The  energy  equation  shows 
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the  balance  between  the  time  rate  of  change  of  energy  and  the  combination  of  com- 
pressional  heating,  divergence  of  heat  flow,  local  heating  and  cooling,  and  transfer  of 
energy  due  to  elastic  collisions. 

Solving  the  coupled,  nonlinear,  continuity,  momentum,  and  energy  equations  for 
each  species  is  required  to  obtain  a  representative  profile  of  the  ion  and  electron  con¬ 
centrations  in  the  ionosphere  and  plasmasphere.  A  numerical  solution  is  needed  as 
this  is  not  possible  analytically  without  making  trivializing  assumptions.  Figure  8 
shows  a  representative  vertical  profile  of  the  major  ions  in  the  ionosphere  and  plas¬ 
masphere  created  using  the  numerical  solution  of  the  IPM.  Molecular  ions  dominate 
in  the  E  region  while  0+  is  the  principal  ion  in  the  F  region  and  topside  ionosphere. 
The  transition  to  the  plasmasphere  occurs  where  H+  becomes  the  primary  ion. 


Log  Electron  Density  Log  Electron  Density 


Figure  8.  Ion  densities  in  the  ionosphere  and  plasmasphere.  Molecular  ions  dominate 
in  the  E  region  while  0+  is  the  principal  ion  in  the  F  region  and  topside  ionosphere.  II 
becomes  the  primary  ion  in  the  plasmasphere  (Adapted  from  Scherliess  et  al.  (2004)) 
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2.3  Ionosphere-Plasmasphere  Model 


The  Ionosphere-Plasmasphere  Model,  developed  at  USU  by  Dr.  Robert  Schunk, 
is  a  physics-based  model  that  calculates  the  three-dimensional  time  dependent  ion 
and  electron  density  distributions.  The  IPM  solves  the  continuity  and  momentum 
equations  for  six  ions  ( NO+ ,  0^~,  0+ ,  H+ ,  and  He+)  on  flux  tubes  that  follow 

the  geomagnetic  field  and  allow  for  inter-hemispheric  flow  ( Schunk  et  al ,  2004).  A 
flux  tube  is  a  cylindrical  region  of  space  around  a  magnetic  field  line.  The  cross- 
sectional  area  of  the  tube  and  the  strength  of  the  magnetic  field  within  the  tube 
may  change,  but  the  magnetic  flux  through  the  tube  always  remains  constant.  The 
use  of  flux  tubes  simplifies  the  calculations  by  allowing  the  plasma  densities  of  one 
flux  tube  to  be  solved  independently  of  the  other  flux  tubes.  Due  to  the  complex 
chemical  reactions,  the  D  region  is  not  included  in  the  IPM.  In  the  E  region,  chemical 
equilibrium  is  assumed  for  NO+ ,  (A%  ,  and  0+ ,  and  the  four  continuity  equa¬ 

tions  are  solved  simultaneously  at  each  grid  point.  Above  the  E  region,  both  the 
continuity  and  momentum  equations  are  solved  for  H+ ,  0+ ,  and  He+.  Included  in 
the  calculations  are  the  effects  of  the  equatorial  electric  field  and  interactions  with 
the  neutral  atmosphere.  Although  solving  along  individual  flux  tubes  creates  output 
that  is  two-dimensional,  calculating  and  tracking  the  ion  densities  for  many  flux  tubes 
across  many  longitudes  results  in  a  three-dimensional  global  model. 

The  IPM  can  be  run  globally,  regionally,  or  locally,  and  covers  geomagnetic  lati¬ 
tudes  from  60°  N  to  60°  S  while  extending  vertically  from  90-30000 km  to  cover  the  E 
and  F  regions  as  well  as  the  topside  ionosphere  and  the  plasmasphere  ( Schunk  et  al, 
2004).  Thompson  et  al.  (2009)  shows  that  inclusion  of  the  plasmasphere  in  the  IPM  is 
important  due  to  model  assimilation  of  GPS  total  electron  content  (TEC)  data.  The 
GPS  path  passes  through  the  entire  plasmasphere,  and  the  plasmasphere’s  nighttime 
contribution  to  the  TEC  can  exceed  50%  ( Scherliess  et  al,  2009). 
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Although  the  IPM  output  is  on  a  geomagnetic  grid,  a  multi-step  transformation 
of  the  field-aligned  equations  is  performed  within  the  model  in  order  to  allow  a  more 
efficient  numerical  solution  ( Schunk  et  al,  2004;  Schunk  and  Nagy,  2009).  First  the 
transport  equations  are  transformed  to  spherical  coordinates,  then  to  dipolar  coordi¬ 
nates,  and  finally  to  a  “sinh”  variable.  This  transformation  optimizes  the  number  of 
grid  points  along  a  field  line.  The  result  is  unequally  spaced  grid  points  along  each 
flux  tube  (Figure  9).  The  number  of  flux  tubes  varies  with  altitude;  from  90-600/cm 
the  flux  tubes  cross  the  equator  at  an  equally  spaced  distance  of  20km.  This  spacing 
increases  exponentially  with  increasing  altitude;  a  similar  variation  in  vertical  resolu¬ 
tion  has  been  successfully  used  in  the  IFM  (Appendix  A).  The  flux  tube  spacing  has 
also  been  optimized  in  order  for  the  model  to  converge  at  the  transition  heights.  As 
vertical  drifts  lift  and  lower  the  ionosphere,  flux  tubes  are  added  to  or  deleted  from 
the  upper  and  lower  boundaries  as  required  to  maintain  altitude  boundaries. 

The  combination  of  unequally  spaced  grid  points  and  flux  tubes  results  in  more 
grid  points  in  the  lower  altitudes  where  more  physical  interactions  occur  and  fewer 
grid  points  at  higher  altitudes  where  densities  are  lower.  Figure  9  shows  the  latitude- 
altitude  grid  used  in  the  IPM  for  both  the  ionosphere  and  plasmasphere.  Latitudinal 
resolution  is  greatest  over  the  equatorial  anomalies  and  around  2.5°  at  mid-latitudes. 

The  longitudinal  resolution  depends  on  the  number  of  longitudinal  planes  solved 
in  the  model.  In  the  global  mode,  a  10°  resolution  is  used  with  36  planes.  The  reso¬ 
lution  can  be  increased  by  decreasing  the  grid  spacing  to  1°  in  the  regional  and  local 
modes;  however,  increasing  the  number  of  planes  significantly  increases  computational 
runtime.  Figure  10  shows  the  horizontal  spacing  of  the  IPM  output  using  48  longi¬ 
tudinal  planes  (7.5°  resolution).  Note  that  in  the  IPM  the  output  is  in  geomagnetic 
coordinates  but  is  plotted  in  geographic  coordinates.  The  International  Geomagnetic 
Reference  Field  (IGRF)  is  used  in  the  IPM  to  model  the  earth’s  magnetic  field. 
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Figure  9.  Latitude-altitude  grid  used  in  the  IPM  showing  the  vertical  and  latitudinal 
resolution  of  the  IPM  for  the  ionosphere  (left)  and  the  plasmasphere  (right)  (Adapted 
from  Scherliess  et  al.  (2004)) 
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Figure  10.  Example  of  longitudinal  resolution  in  the  IPM  using  48  longitudinal  planes. 
Plasma  densities  at  300 km  at  1200UT  during  solar  medium  conditions  are  displayed. 
The  ‘S’  shape  of  the  longitudinal  planes  is  due  to  the  geomagnetic  flux  tubes  graphed 
geographically  (Adapted  from  Scherliess  et  al.  (2004)) 
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The  IPM  calculates  the  ionospheric  drivers  needed  to  run  the  model  using  internal 
empirical  models.  The  Mass  Spectrometer  and  Incoherent  Scatter  (MSIS)  model  is 
used  for  the  neutral  density  and  temperature  profiles.  This  empirical  model  provides 
these  parameters  based  on  a  given  set  of  geophysical  conditions  (solar  cycle,  season, 
geomagnetic  activity,  etc.)  ( Hedin ,  1991).  The  zonal  and  meridional  neutral  winds  are 
obtained  from  the  Horizontal  Wind  Model  (HWM)  developed  by  Hedin  et  al.  (1991). 
The  dynamo  electric  fields  are  calculated  in  the  IPM  using  the  model  developed 
by  Fejer  et  al.  (1999).  The  Titheridge  (1998)  model  is  used  in  the  IPM  to  solve  the 
energy  equation  and  calculate  the  ion  and  electron  temperature  profiles. 

Each  empirical  model  is  solved  independently  within  the  IPM;  in  addition,  the 
E  region,  Fi  region,  and  F2  region  are  also  solved  in  separate  modules  within  the 
IPM.  As  a  final  step,  the  IPM  splices  together  the  solutions  for  the  different  regions. 
Because  each  module  is  run  separately  within  the  IPM,  the  solutions  may  not  be  self- 
consistent.  For  example,  although  the  zonal  thermospheric  winds  play  an  important 
role  in  setting  up  the  equatorial  electric  fields,  the  winds  and  the  electric  fields  are 
calculated  separately  with  their  respective  empirical  models;  therefore,  the  HWM 
calculated  zonal  wind  does  not  affect  the  electric  fields  ( Schunk ,  2010).  The  zonal 
wind  does,  however,  play  a  role  in  transport  of  the  plasma  because  the  momentum 
equation  requires  this  term  as  part  of  its  solution. 

Using  empirical  model  results  to  solve  the  physical  equations  in  the  physics-based 
model  can  lead  to  erroneous  and  inconsistent  features  in  the  physics-based  model 
output.  Data  assimilation  can  be  used  to  correct  these  errors,  however,  if  data  as¬ 
similation  is  not  accomplished  or  if  there  is  only  a  limited  amount  of  data  the  model 
must  be  corrected  by  comparing  the  output  to  a  large  set  of  observational  data  and 
manually  adjusting  physical  parameters  in  the  model  prior  to  its  use.  The  result  of 
correcting  these  errors  is  a  more  realistic  plasma  density  distribution. 
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III.  Methodology 


3.1  Chapter  Overview 

This  chapter  describes  the  methods  used  to  perform  this  study  and  is  broken  into 
five  sections.  The  first  section  describes  the  interpolation  of  the  IPM  from  magnetic 
flux  tubes  to  a  global  geographic  grid.  The  selection  of  the  geophysical  conditions 
chosen  for  this  study  is  discussed  in  section  two.  The  third  section  discusses  the 
model  runs  including  details  on  the  model  inputs  and  the  adjustments  made  for  each 
comparison.  The  fourth  section  outlines  how  the  variables  were  calculated  from  the 
IPM  output  and  plotted.  The  final  section  describes  the  techniques  used  to  compare 
the  model  output. 

3.2  IPM  Interpolation 

Before  the  model  physics  could  be  investigated,  the  IPM  had  to  be  interpolated 
to  output  on  a  three-dimensional  global  grid.  The  IPM  interpolation  code  takes 
the  IPM  output  from  a  two-dimensional  longitude  versus  altitude  geomagnetic  grid 
consisting  of  many  flux  tubes  (Figure  9)  and  converts  it  to  a  global  three-dimensional 
(latitude,  longitude,  altitude)  geographic  grid.  The  grid  chosen  was  a  3.5°  latitude 
by  7.5°  longitude  grid  with  variable  altitude  spacings  (Figure  11).  The  interpolation 
code  searches  through  each  longitudinal  plane  within  a  set  number  of  degrees  of  the 
desired  latitude  grid  point  and  fills  a  vertical  array  with  values  going  from  low  to  high 
altitude.  This  array  is  then  fit  with  a  cubic  spline  interpolation  routine.  Once  all  of 
the  latitudes  have  been  filled,  a  cubic  spline  in  longitude  is  performed  to  match  the 
data  to  the  longitudes  of  the  IFM  output  grid  ( Gardner ,  2010).  The  longitudes  have 
to  be  fit  as  well  because  the  IPM  flux  tubes  form  a  type  of  ‘S’  shape  in  geographic 
coordinates  (see  Figure  10). 
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Figure  11.  Vertical  resolution  of  the  interpolated  IPM  output  (not  to  scale) 


Multiple  comparisons  of  the  interpolated  IPM  output  to  the  IFM  output  ensured 
that  the  IPM  was  producing  valid  results  and  that  the  interpolation  code  was  working 
correctly.  A  general  agreement  between  the  IFM  and  IPM  was  achieved,  however, 
the  model  output  differed  in  some  instances  due  to  the  differences  in  each  model’s 
physics.  The  physical  differences  between  the  IFM  and  IPM  are  given  in  Table  1.  The 
interpolation  code  for  the  IPM  was  written  by  USU  and  modified  with  the  aid  of  the 
IFM/IPM  comparisons  completed  at  the  Air  Force  Institute  of  Technology  (AFIT) 
for  this  thesis. 


Table  1.  Differences  in  the  IFM  and  IPM 


IFM 

IPM 

Simple  magnetic  field 

IGRF  magnetic  field 

Simple  H+  calculation 

Rigorous  H+  calculation 

No  He+ 

He+  included 

Coarse  geographic  grid 

Geomagnetic  flux  tube  grid 

1600fcm  top  boundary 

30000 km  top  boundary 

Temperatures  numerically  solved 

Empirically-derived  temperatures 
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3.3  Geophysical  Condition  Selection 


The  ionosphere  varies  significantly  with  season  and  solar  cycle.  In  order  to  cap¬ 
ture  the  large  variation  in  these  conditions  while  limiting  the  number  of  times  the 
model  needed  to  be  run,  three  sets  of  geophysical  conditions  were  chosen:  solar  maxi¬ 
mum  June  solstice,  solar  minimum  December  solstice,  and  solar  medium  fall  equinox. 
The  combination  of  solar  maximum  and  June  solstice  results  in  the  highest  levels  of 
ionization  in  the  northern  hemisphere  due  to  the  high  solar  activity  and  increased 
daylight  during  summer.  High  solar  activity  during  solar  maximum  corresponds  to 
increased  solar  EUV  radiation  which  is  absorbed  in  earth’s  ionosphere  through  in¬ 
creased  ionization.  Solar  minimum  December  solstice  conditions  have  the  opposite 
result:  decreased  solar  activity  with  less  solar  EUV  radiation  and  decreased  daylight. 
Solar  medium  fall  equinox  was  also  chosen  as  a  middle  ground  with  medium  levels 
of  solar  activity  and  equal  daylight  in  each  hemisphere.  Table  2  summarizes  the 
geophysical  conditions  for  which  the  models  were  run. 


Table  2.  Geophysical  conditions  chosen  for  the  model  runs 


Case  # 

Solar  Cycle 

Season 

Year 

Day 

-Fio.7 

KP 

Ap 

1 

Minimum 

December  Solstice 

1995 

358 

70 

3 

15 

2 

Medium 

Fall  Equinox 

1998 

267 

140 

3 

15 

3 

Maximum 

June  Solstice 

2001 

172 

220 

3 

15 

3.4  Model  Runs 

The  IPM  requires  as  input  the  day,  year,  Universal  Time  (UT)  start  time,  output 
interval,  daily  and  90-day  average  F10  7  solar  flux,  and  daily  Ap  index.  In  addition,  the 
resolution  and  domain  of  the  IPM  are  set  by  specifying  the  number  of  degrees  between 
each  longitudinal  plane  and  the  upper  and  lower  altitude  boundaries.  The  altitudes 
are  in  reference  to  the  height  at  which  the  flux  tubes  cross  the  equatorial  plane. 
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For  this  project  the  IPM  was  set  to  output  36  longitudinal  planes  (10°  longitudinal 
resolution)  with  altitude  bounds  of  135-20000/cm  at  the  equator. 

Once  the  model  run  has  been  started,  four  days  of  output  are  required  before 
the  model  reaches  a  steady  state  solution.  For  example,  if  the  model  is  started  on 
day  354,  the  output  will  not  be  valid  until  day  358.  To  ensure  that  the  model  has 
indeed  converged  to  a  steady  state  solution,  the  fourth  day  of  output  (day  358)  was 
compared  to  the  fifth  day  of  output  (day  359)  (Figure  12,  top  row)  for  the  baseline 
run  of  case  1.  In  addition,  the  fourth  day  was  also  compared  to  the  twentieth  day 
(day  374)  to  show  that  the  model  first  converges  and  then  outputs  the  same  steady 
state  solution  from  that  point  on  (Figure  12,  bottom  row).  For  these  comparisons 
the  input  parameters  were  not  changed.  Although  only  the  0000UT  TEC  is  shown 
in  Figure  12,  the  comparisons  were  done  for  all  24  hours  and  included  variables  such 
as  NmF2,  hmF2,  300 km  Ne,  and  800/cm  Ne.  The  left  hand  plots  of  Figure  12  show 
the  output  from  the  fourth  day  of  the  model  run,  the  middle  plots  show  the  fifth 
day  (top)  and  twentieth  day  (bottom),  and  the  right  hand  plots  show  the  difference. 
As  can  be  seen  in  both  difference  plots,  there  is  less  than  a  ±1  TECU  change  in 
either  of  the  comparisons.  The  cause  for  the  differences  in  the  bottom  right  plot 
around  60°E  and  300°E  are  not  known  but  are  not  in  regions  that  are  evaluated  in 
this  project  and  therefore  do  not  affect  the  results.  Further  investigation  of  the  cause 
for  the  differences  in  these  regions  is  recommended  for  future  work.  The  negligible 
differences  across  the  rest  of  the  plot,  which  can  be  attributed  to  solar  zenith  angle 
changes,  confirms  the  assumption  that  the  model  has  converged  to  the  steady  state 
solution  in  the  areas  of  interest  by  the  fourth  day  of  output. 

With  steady  state  conditions  not  reached  until  the  fourth  day  of  output  and  the 
1:1  ratio  of  model  runtime  to  output,  the  IPM  takes  at  least  five  days  of  runtime 
to  obtain  one  day  of  usable  output.  To  force  the  model  to  run  faster,  the  IPM  was 
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t  igure  12.  1EC  steady  state  comparison  tor  the  tour  day  run  compared  to  a  hve  day  run  (top)  and  the  tour  day  run  compared 
to  a  twenty  day  run  (bottom).  The  left  plot  is  the  four  day  run,  the  middle  plot  is  the  5  and  20  day  run,  and  the  right  plot  is 
the  difference.  After  day  five,  the  model  output  has  reached  steady  state  conditions  and  differs  by  less  than  1  TECU,  which 
can  be  attributed  to  the  change  in  zenith  angle 


ran  with  three  different  altitude  sections  (135-5000fcm,  5000-1 1000km,  and  11000- 
20000A’m) .  These  sections  are  recombined  into  one  continuous  altitude  section  using 
the  interpolation  code.  The  IPM  does  not  require  an  end  time  and  will  continue 
until  the  run  is  manually  terminated.  The  IPM  outputs  the  electron  density,  six  ion 
densities,  and  temperature  arrays  in  large  text  files.  The  interpolation  code  can  be 
run  on  any  of  these  output  variables. 

To  investigate  the  physical  parameters  in  the  IPM,  an  input  file  was  created  to 
turn  ‘off’  and  ‘on’  different  sections  of  code  in  the  IPM.  Turning  ‘on’  a  section  of 
code  changed  the  parameter  specified  in  the  input  file.  The  adjustments  for  each  of 
the  five  parameters  ( 0+/0  collision  frequency,  daytime  production  rate,  nighttime 
drift,  zonal  winds,  and  tidal  forcing),  in  addition  to  the  default  settings,  are  given  in 
Table  3  and  are  discussed  in  more  detail  in  Chapter  IV.  These  physical  parameters  are 
known  to  be  uncertain  in  physics-based  models  and  the  adjustments  were  motivated 
by  comparisons  with  the  12-year  TOPEX  TEC  database  ( Schunk ,  2010). 


Table  3.  Adjustments  made  to  the  IPM 


Parameter 

Default 

Adjustment 

0+/0  Collision 
Frequency 

Normal 

Doubled 

Zonal  Winds 

HWM  derived  zonal  winds 

Zonal  winds  set  to  zero 

Tidal  Forcing 

No  tidal  forcing 

Tidal  forcing  included  by 
modulating  the  E  x  B 
drift 

Daytime 

Production 

Production  multiplied  by  1.8 

Multiplication  factor  de¬ 
creased  as  a  linear  func¬ 
tion  of  Fio.7 

Nighttime  Drifts 

No  correction  to  drifts 

Decreased  the  nighttime 
E  x  B  drift  as  a  linear 
function  of  F10.7 

In  total,  18  model  runs  were  completed:  1  baseline  run  and  5  adjustments  for 
3  geophysical  cases.  The  model  runs  were  done  in  collaboration  with  Utah  State 


27 


University.  USU  developed  both  the  IFM  and  IPM  and  changes  to  the  models  were 
made  by  USU  researchers.  The  models  were  run  remotely  on  USU  computers;  the 
model  output  was  downloaded  to  AFIT  computers  where  the  analysis  was  completed. 

3.5  Variable  Calculation  and  Data  Visualization 

Data  visualization  is  an  important  part  of  research  because  it  allows  large  arrays 
of  data  to  be  interpreted  quickly  and  easily.  This  project  produced  hundreds  of  arrays 
of  data  while  investigating  the  total  electron  content  (TEC),  electron  density  (Ne), 
and  F2  region  peak  density  (iVmF2)  and  height  (hm.F2)  of  the  ionosphere.  This  section 
discusses  how  these  variables  were  calculated  and  plotted  from  the  IPM  output,  and 
the  technique  used  to  create  the  longitude  slice  plots  and  the  local  time  plots. 

The  IPM  outputs  a  3-D  (latitude,  longitude,  altitude)  array  of  electron  densities 
for  every  hour  (in  UT)  in  a  text  hie.  This  data  was  read  into  MATLAB  where  it  was 
then  manipulated  to  compute  the  required  variables  and  to  plot  the  data  in  a  variety 
of  configurations.  The  TEC  was  calculated  by  adding  up  the  electron  density  at  every 
altitude  for  each  latitude/longitude  pair  and  then  converting  the  densities  to  TECU 
by  dividing  by  1016  electrons/m2.  The  NmF2  and  hmF2  were  calculated  together 
by  systematically  stepping  through  every  altitude  and  locating  the  largest  value  for 
every  latitude/longitude  pair  (NmF2)  and  the  height  at  which  that  value  occurred 
(hmF2).  These  variables  were  saved  in  2-D  arrays  for  every  hour  calculated.  The 
electron  density  arrays  at  300 km,  400 km,  and  800/cm  were  also  saved  in  2-D  arrays 
for  comparison  and  evaluation.  The  longitude  slice  plots  were  created  in  a  manner 
similar  to  the  Ne  plots.  Instead  of  choosing  one  specific  altitude,  one  longitude  was 
selected  to  display.  The  longitude  slice  plots  are  more  intuitive  when  plotted  in  local 
time;  this  was  accomplished  by  simply  calculating  the  local  time  of  each  plot  by 
utilizing  the  UT  offset  for  the  chosen  longitude. 
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The  2-D  arrays  were  plotted  using  MATLAB’s  imagesc  plotting  function  and  a 
72-hue  colormap.  The  minimum  value  and  maximum  value  of  the  entire  24-hour 
period  were  found  for  each  variable  and  these  values  were  used  as  the  boundaries  of 
the  colorscale.  This  ensured  that  the  colorscale  did  not  vary  between  plots.  The  IPM 
does  not  cover  the  high  latitude  regions  and  so  these  areas  (poleward  of  ±60°)  were 
left  white  to  indicate  no  data. 

Although  the  UT  plots  are  the  most  intuitive  (day  versus  night  is  easy  to  distin¬ 
guish  and  understand) ,  plotting  in  local  time  can  highlight  diurnal  trends  in  the  data. 
When  the  same  local  time  for  each  location  is  plotted  together,  the  variations  due 
to  geographic  location  (instead  of  time)  can  easily  be  seen.  For  example,  when  the 
TEC  is  plotted  at  1200LT  (Figure  13),  the  longitudinal  variation  in  the  equatorial 
anomalies  is  easily  identifiable.  This  type  of  plot  was  created  by  determining  the  UT 
offset  for  each  longitude  and  then  choosing  the  correct  local  time  longitude  from  each 
UT  array  and  saving  the  data  to  a  new  local  time  array.  This  array  was  then  plotted 
in  the  same  manner  as  the  UT  arrays. 


Total  Electron  Content 
Solar  Med  Fall  Equinox  (1998  267)  1200LT 


TEC  Units 

Figure  13.  TEC  plotted  in  local  time  highlights  the  variation  in  TEC  across  longitudes 
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3.6  Analysis  Techniques 


The  next  step  in  this  project  was  to  evaluate  how  adjusting  the  physical  parame¬ 
ters  affected  the  model  output.  For  this  evaluation,  two  comparison  techniques  were 
used — a  difference  and  a  percentage  change.  These  comparisons  were  done  over  the 
entire  spatial  grid  in  order  to  highlight  regions  of  differences.  The  difference  was 
computed  by  subtracting  the  default  IPM  output  value  at  each  grid  point  from  the 
adjusted  IPM  output  value  at  the  corresponding  grid  point  (Equation  17).  The  per¬ 
cent  increase  (or  decrease)  was  calculated  by  dividing  the  difference  by  the  default 
IPM  output  value  and  then  multiplying  by  100  to  get  a  percentage  (Equation  18). 

Difference  IPMac^US(ec/  IPMt/eja(i^  (17) 

t,  i  T  IPM  adjusted,  IP  default  i  r>n  /-.  0\ 

Percent  Increase  =  - TT^  - - -  x  100  (18) 

IP  Mde  fault  K  ' 

These  comparisons  were  done  for  each  variable  of  interest  (TEC,  NmF2,  hmF2, 
and  400 km  Ne)  every  hour  for  24  hours  and  plotted  in  both  UT  and  local  time 
(LT).  Additionally,  different  longitudes  were  chosen  to  do  a  latitude  versus  altitude 
(longitude  slice)  comparison  of  the  electron  density.  In  general,  the  longitudes  were 
selected  based  on  where  the  TEC  or  NmF2  showed  large  changes.  The  longitude  slices 
are  displayed  only  in  local  time. 
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IV.  Results  and  Analysis 


4.1  Chapter  Overview 

The  objective  of  this  project  was  to  complete  a  sensitivity  analysis  of  certain 
physical  parameters  in  a  physics-based  model  and  to  examine  how  adjusting  these 
parameters  affects  the  model  output.  The  parameters  that  were  investigated  are 
known  to  have  inherent  uncertainty  within  models.  Each  of  the  parameters  were 
adjusted  independently  of  each  other  to  isolate  the  effects  of  each  adjustment  on  the 
model  output.  First  the  baseline  results  for  each  geophysical  case  (see  Chapter  II)  are 
given  where  no  adjustments  were  made  to  the  model.  Next,  the  effects  of  adjusting 
each  parameter  are  discussed.  The  0+/0  collision  frequency  is  examined  first  followed 
by  the  zonal  winds,  daytime  production,  nighttime  E  x  B  drift,  and  finally  the  four- 
wave  tidal  influence.  Comparing  the  baseline  output  to  the  adjusted  output  gives  an 
understanding  of  the  importance  of  correctly  specifying  these  physical  parameters  in 
a  physics-based  model. 

4.2  Baseline  Runs 

This  section  gives  the  default  output  from  the  IPM  for  which  no  adjustments 
were  made  to  any  parameters.  Understanding  the  baseline  results  is  imperative  to 
understanding  how  adjusting  the  physical  parameters  affects  this  output.  Figure  14 
shows  plots  of  TEC  for  all  three  geophysical  conditions  at  both  0000UT  and  1200UT. 
Maximum  daytime  values  of  TEC  range  from  45-160  TECU  while  nighttime  values 
are  as  low  as  2-3  TECU.  As  expected,  TEC  is  a  function  of  the  incoming  solar  flux  and 
increases  from  solar  minimum  to  solar  maximum.  In  the  0000UT  plots,  1200LT  is  near 
180°U  and  the  higher  daytime  TEC  values  are  due  to  the  presence  of  photoionization. 
The  terminators  can  also  be  clearly  seen  in  these  plots  with  the  area  of  illumination 


31 


Total  Electron  Content  for  All  Three  Geophysical  Cases 
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Figure  14.  TEC  for  all  three  geophysical  cases  (solar  minimum  December  solstice  - 
top,  solar  medium  fall  equinox  -  middle,  and  solar  maximum  June  solstice  -  bottom) 
at  both  OOOOUT  (left)  and  1200UT  (right).  Both  solar  cycle  and  seasonal  effects  can 
be  seen  in  the  differences  in  TEC  in  these  plots.  The  circled  region  is  the  Weddell  Sea 
Anomaly.  Note  that  the  scales  are  different  for  each  geophysical  case 
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changing  according  to  season.  The  equatorial  anomaly  regions  of  enhanced  TEC  on 
each  side  of  the  geomagnetic  equator  are  also  evident  in  these  plots. 

Similar  to  the  TEC  plots,  the  NmF2  for  all  three  geophysical  conditions  is  plotted 
in  Figure  15.  These  plots  show  the  peak  density  in  the  F2  region  of  the  ionosphere, 
which  is  a  significant  contributor  to  the  overall  TEC.  Because  of  this,  the  same  features 
(equatorial  anomalies,  and  day/night,  seasonal,  and  solar  cycle  differences)  can  also 
be  seen  in  these  plots.  The  scale  of  these  plots  is  the  log  of  the  electron  density,  i.e. 
a  value  of  ‘5’  means  105  electrons/cm-3. 

Also  of  note  in  the  baseline  runs  is  the  area  of  enhanced  TEC  and  NmF2  over  the 
Southeast  Pacific  during  solar  minimum  and  solar  medium  (circled  in  the  top  and 
middle  left  hand  plots  of  Figures  14  and  15).  This  feature  shows  up  throughout  the 
comparisons  and  is  a  result  of  the  Weddell  Sea  Anomaly  (WSA).  The  WSA  in  the 
ionosphere  is  characterized  by  higher  plasma  densities  at  night  than  during  the  day 
in  the  region  near  the  Weddell  Sea  (Antarctica).  This  nighttime  enhancement  is  due 
to  the  vertical  plasma  drifts  resulting  from  the  combination  of  thermospheric  neutral 
winds  and  the  magnetic  held  configuration  in  this  region  ( Horvath  and  Essex,  2003). 
Jee  et  al.  (2009)  found  that  the  WSA  is  most  prominent  during  December  solstice 
(southern  hemisphere  summer)  and  occurs  in  all  seasons  for  low  FW  7  values  while 
for  high  F10.7  values  it  is  not  as  prominent  during  southern  hemisphere  winter  (June 
solstice) . 

The  altitude  structure  of  the  electron  density  is  also  important  because  densities 
can  increase  at  one  altitude  while  decreasing  at  another  with  no  resulting  change  in 
TEC.  Figure  16  shows  a  latitude  versus  altitude  plot  of  the  the  electron  density  at 
90°  E  at  both  0000LT  and  1200LT.  Both  the  solar  cycle  and  seasonal  dependence 
of  Ne  is  evident  in  these  plots.  As  solar  activity  increases,  short  wavelength  solar 
radiation  increases  resulting  in  higher  levels  of  photoionization  and  Ne.  The  seasonal 
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F2  Peak  Electron  Density  for  All  Three  Geophysical  Cases 
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Figure  15.  NrnF2  for  all  three  geophysical  cases  (solar  minimum  December  solstice  - 
top,  solar  medium  fall  equinox  -  middle,  and  solar  maximum  June  solstice  -  bottom) 
at  both  OOOOUT  (left)  and  1200UT  (right).  Both  solar  cycle  and  seasonal  effects  can 
be  seen  in  the  differences  in  NrnF2  in  these  plots.  The  circled  region  is  the  Weddell  Sea 
Anomaly.  Note  that  the  scales  are  different  for  each  geophysical  case 
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90°E  Electron  Density  for  All  Three  Geophysical  Cases 
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Figure  16.  Ne  at  90°E  for  all  three  geophysical  cases  (solar  minimum  December  solstice 
-  top,  solar  medium  fall  equinox  -  middle,  and  solar  maximum  June  solstice  -  bottom) 
at  both  0000UT  (left)  and  1200UT  (right).  Both  solar  cycle  and  seasonal  effects  can 
be  seen  in  the  differences  in  Ne  in  these  plots.  Note  that  the  scales  are  different  for 
each  geophysical  case 
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dependence  is  seem  in  the  winter  anomaly — higher  plasma  densities  in  the  winter 
hemisphere  than  in  the  summer  hemisphere.  Even  though  the  ion  production  rate  is 
higher  in  the  summer  hemisphere,  the  loss  rate  is  also  much  higher  and  overwhelms 
the  increase  in  ionization  due  to  production.  Additionally,  the  plasma  is  transported 
by  meridional  neutral  winds  along  the  held  lines  from  the  summer  (hot)  hemisphere 
to  the  winter  (cold)  hemisphere.  These  processes  result  in  the  summer  hemisphere 
density  peak  being  at  a  higher  altitude  and  lower  in  magnitude  than  the  winter 
hemisphere  density  peak.  This  structure  is  evident  in  the  top  and  bottom  plots  of 
Figure  16.  These  plots  are  especially  useful  to  show  the  development  and  subsequent 
depletion  of  the  equatorial  anomalies  as  well  as  the  location  and  magnitude  of  the 
peak  density. 

4.3  0+/0  Collision  Frequency 

The  0+/0  collision  frequency  was  the  first  parameter  investigated.  For  this  com¬ 
parison,  the  standard  value  of  the  0+/0  collision  frequency  as  calculated  by  Schunk 
and  Nagy  (2009)  was  doubled.  This  directly  affects  both  diffusion  along  the  magnetic 
field  and  the  plasma’s  interaction  with  the  thermospheric  neutral  wind  (see  Figure  5). 
The  diffusion  rate  is  inversely  proportional  to  the  collision  frequency,  therefore  dou¬ 
bling  the  collision  frequency  reduces  diffusion  by  a  factor  of  two.  The  result  is  slower 
downward  diffusion  of  0+  which  aids  in  maintaining  the  F  region  plasma.  Addition¬ 
ally,  doubling  the  0+/0  collision  frequency  acts  to  significantly  increase  the  coupling 
of  the  plasma  drifts  to  the  neutral  wind.  The  thermospheric  neutral  winds  are  di¬ 
rected  radially  away  from  the  subsolar  point  with  daytime  meridional  winds  blowing 
toward  the  poles  and  nighttime  meridional  winds  blowing  toward  the  equator.  During 
the  day  the  wind  becomes  more  effective  in  driving  the  plasma  to  lower  altitudes,  and 
at  night  the  wind  is  more  effective  in  driving  the  plasma  to  higher  altitudes. 
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During  the  day,  the  reduced  0+  diffusion  and  the  increased  neutral  wind  effect  act 
to  oppose  each  other.  Even  though  the  neutral  winds  are  doubled,  daytime  neutral 
winds  are  weak,  and  so  the  wind  effect  and  diffusion  effect  are  nearly  equal  resulting  in 
very  little  change  in  the  peak  electron  density.  However,  if  the  wind  effect  is  stronger 
than  the  diffusion  effect,  a  decrease  in  peak  electron  density  may  be  observed. 

At  night,  the  diffusion  and  the  winds  act  together  to  maintain  F  region  plasma 
densities.  The  overall  effect  is  higher  peak  densities  at  higher  altitudes  because  the 
plasma  is  lifted  to  higher  altitudes  where  recombination  rates  are  slower.  During 
sunrise  and  sunset,  the  neutral  winds  are  primarily  in  the  east- west  direction  across 
the  terminator  and  do  not  significantly  affect  the  plasma  densities,  except  in  regions 
where  the  magnetic  declination  is  large. 

Figure  17  shows  the  percent  increase  in  NmF2  when  the  0+/0  collision  frequency 
is  doubled.  The  0000UT  (left  hand  plots)  and  1200UT  (right  hand  plots)  are  shown 
for  all  three  geophysical  cases.  The  same  general  changes  occur  for  all  three  cases 
with  a  significant  increase  in  NmF2  at  night  and  very  little  change  in  NmF2  during 
the  day.  The  largest  increase  is  250%  and  occurs  in  the  equatorial  anomalies  during 
solar  minimum  December  solstice  (top  plots)  due  to  the  combination  of  the  HWM 
derived  neutral  winds  and  the  configuration  of  the  magnetic  field.  In  regions  where 
the  neutral  wind  efffect  dominates  over  the  diffusion  effect  during  the  day,  up  to  a 
30%  decrease  can  be  seen  in  the  NmF2. 

In  order  to  investigate  more  closely  what  is  occurring  when  the  0+/0  collision 
frequency  is  doubled,  mid-latitude  locations  at  four  equally  spaced  longitudes  in  both 
the  northern  and  southern  hemisphere  were  chosen  and  the  peak  electron  density 
for  the  normal  0+/0  collision  frequency  as  well  as  for  the  doubled  0+/0  collision 
frequency  were  plotted  versus  local  time  (Figures  18-21).  The  expected  trend  appears 
in  all  plots  with  the  greatest  increase  in  NmF2  occurring  at  night  and  a  minimal 


37 


NmF2  Percent  Increase  with  a  Doubled  0+/0  Collision  Frequency 
for  All  Three  Geophysical  Cases 
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Figure  17.  NrnF2  increase  with  doubled  0+/0  collision  frequency  for  all  three  geophysical 
cases  (solar  minimum  December  solstice  -  top,  solar  medium  fall  equinox  -  middle,  and 
solar  maximum  June  solstice  -  bottom)  at  both  0000UT  (left)  and  1200UT  (right).  The 
largest  changes  in  NrnF2  occur  during  solar  minimum.  Note  that  the  scales  are  different 
for  each  geophysical  case 


increase  occurring  during  the  day.  It  was  found  that  the  peak  increase  occurs  in  the 
early  morning  hours  after  plasma  densities  have  had  time  to  build  up  at  the  higher 
altitudes.  The  plasma  decay  at  sunset  was  found  to  be  slower  with  the  doubled 
0+/0  collision  frequency  due  to  the  reduced  downward  diffusion.  The  same  general 
behavior  is  seen  for  all  eight  locations  with  the  largest  changes  occurring  during  solar 
minimum  December  solstice  and  solar  medium  fall  equinox  at  45°  S  90°  E  and  45°  S 
180°  E.  In  these  instances,  the  NmF2  increases  by  120-180%.  The  nonlinearity  of  the 
changes  is  a  result  of  the  complex  nature  of  the  many  interactions  in  the  ionosphere. 
Nightime  peak  heights  increased  by  10-50/cm  while  daytime  heights  only  increased  a 
maximum  of  10 km. 

Of  interest  are  the  double  peaks,  indicated  by  black  arrows,  in  the  270° E  southern 
hemisphere  NmF2  plots  in  Figure  21.  The  peaks  occur  just  after  sunrise  and  before 
sunset.  This  double  peak  in  NmF2  is  a  result  of  the  Weddell  Sea  Anomaly  discussed  in 
Section  4.2.  The  increased  efficiency  of  the  neutral  winds  during  the  day  act  to  bring 
the  plasma  down  the  field  lines  to  regions  of  higher  recombination  and  thus  lowering 
the  plasma  density.  At  night  the  opposite  occurs  and  the  winds  lift  the  plasma  up  the 
field  lines  resulting  in  enhanced  densities  ( Horvath  and  Essex ,  2003).  As  mentioned 
previously,  the  peaks  are  most  prominent  during  solar  minimum  December  solstice 
and  solar  medium  fall  equinox. 

Table  4  gives  values  for  the  percent  change  in  A%F2  at  local  midnight,  before 
sunrise,  at  local  noon,  and  during  the  daytime  for  all  eight  locations.  As  mentioned 
previously,  the  largest  changes  in  NmF2  occur  at  night,  specifically  just  before  sunrise, 
and  minimal  changes  occur  during  the  daytime.  The  negative  daytime  values  that 
occur  at  solar  maximum  are  a  result  of  a  strong  interaction  with  the  neutral  wind 
that  overwhelms  the  reduced  diffusion  and  drives  the  plasma  down  the  field  lines  to 
where  recombination  is  faster. 
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NmF2  vs.  Local  Time  with  a  Doubled  0+/0  Collision  Frequency 
for  All  Three  Geophysical  Cases 

45°N  0°E  Solar  Min  Dec  Solstice  (1995  358)  45°S  0°E  Solar  Min  Dec  Solstice  (1995  358) 
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Figure  18.  NrnF2  at  0°E  with  doubled  0+/0  collision  frequency  for  all  three  geophysical 
cases  (solar  minimum  December  solstice  -  top,  solar  medium  fall  equinox  -  middle, 
and  solar  maximum  June  solstice  -  bottom)  at  both  45°N  (left)  and  45°S  (right).  The 
black  line  is  the  unchanged  0+/0  collision  frequency,  the  red  line  is  the  doubled  0+/0 
collision  frequency,  and  the  blue  line  is  the  percent  increase.  The  largest  increases 
occur  at  night  with  minimal  changes  during  the  day 
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NmF2  vs.  Local  Time  with  a  Doubled  0+/0  Collision  Frequency 
for  All  Three  Geophysical  Cases 

45°N  90°E  Solar  Min  Dec  Solstice  (1995  358)  45°S  90°E  Solar  Min  Dec  Solstice  (1995  358) 
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Figure  19.  NrnF2  at  90°E  with  doubled  0+/0  collision  frequency  for  all  three  geophysical 
cases  (solar  minimum  December  solstice  -  top,  solar  medium  fall  equinox  -  middle,  and 
solar  maximum  June  solstice  -  bottom)  at  both  45°N  (left)  and  45° S  (right).  The  black 
line  is  the  unchanged  0+/0  collision  frequency,  the  red  line  is  the  doubled  0+/0  collision 
frequency,  and  the  blue  line  is  the  percent  increase.  The  largest  increases  occur  at  night 
with  minimal  changes  during  the  day 


41 


NmF2  vs.  Local  Time  with  a  Doubled  0+/0  Collision  Frequency 
for  All  Three  Geophysical  Cases 
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Figure  20.  NrnF2  at  180°E  with  doubled  0+/0  collision  frequency  for  all  three  geophysical 
cases  (solar  minimum  December  solstice  -  top,  solar  medium  fall  equinox  -  middle,  and 
solar  maximum  June  solstice  -  bottom)  at  both  45°N  (left)  and  45° S  (right).  The  black 
line  is  the  unchanged  0+/0  collision  frequency,  the  red  line  is  the  doubled  0+/0  collision 
frequency,  and  the  blue  line  is  the  percent  increase.  The  largest  increases  occur  at  night 
with  minimal  changes  during  the  day 
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7VmF2  vs.  Local  Time  with  a  Doubled  0+/0  Collision  Frequency 
for  All  Three  Geophysical  Cases 

45°N  270°E  Solar  Min  Dec  Solstice  (1995  358)  45°S  270°E  Solar  Min  Dec  Solstice  (1995  358) 


45°N  270°E  Solar  Med  Fall  Equinox  (1998  267) 
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45°N  270°E  Solar  Max  June  Solstice  (2001  172) 


45°S  270°E  Solar  Max  June  Solstice  (2001  172) 


COFREQ  =  1.0  - COFREQ  =  2.0  - PERCENT  INCREASE 


Figure  21.  NrnF2  at  270°E  with  doubled  0+/0  collision  frequency  for  all  three  geophysical 
cases  (solar  minimum  December  solstice  -  top,  solar  medium  fall  equinox  -  middle,  and 
solar  maximum  June  solstice  -  bottom)  at  both  45°N  (left)  and  45° S  (right).  The  black 
line  is  the  unchanged  0+/0  collision  frequency,  the  red  line  is  the  doubled  0+/0  collision 
frequency,  and  the  blue  line  is  the  percent  increase.  The  largest  increases  occur  at  night 
with  minimal  changes  during  the  day.  Double  peaks  due  to  the  Weddell  Sea  Anomaly 
are  indicated  by  black  arrows 
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Plasma  densities  below  the  peak  were  also  affected  by  doubling  the  0+/ O  collision 
frequency.  The  plasma  in  this  region  was  significantly  depleted  at  sunset  and  filled 
back  in  after  sunrise.  The  upward  motion  of  the  plasma  resulting  from  the  combina¬ 
tion  of  the  increased  wind  effect  and  decreased  0+  diffusion  acts  to  deplete  0+  at  low 
altitudes  with  no  source  of  replenishment.  Figure  22  shows  latitude  versus  altitude 
plots  at  0°E  of  the  percent  increase  in  Ne  due  to  doubling  the  0+/0  collision  fre¬ 
quency  during  solar  minimum  December  solstice  conditions.  The  upper  left  plot  is  at 
midnight  local  time  and  shows  a  nearly  uniform  80%  decrease  in  Ne  from  200-250/.’m. 
As  sunrise  approaches  (upper  middle  and  upper  right  plots),  the  depleted  region  is 
filled  in  as  the  plasma  begins  to  move  back  down  the  field  lines.  By  0800LT  (upper 
right  plot)  the  depletion  is  only  10-20%.  This  slight  depletion  continues  throughout 
the  day  and  then  as  sunset  approaches  and  the  plasma  begins  to  move  to  higher 
altitudes  the  depletion  increases  again.  The  bottom  three  plots  in  Figure  22  show 
the  time  sequence  of  the  density  decrease  during  sunset  from  1800LT-2000LT.  The 
increase  in  densities  above  the  peak  can  also  be  seen  in  these  plots.  These  results  are 
consistent  for  all  three  geophysical  cases. 

4.4  Zonal  Winds 

The  thermospheric  low  latitude  zonal  winds  are  calculated  in  the  IPM  using  the 
HWM.  These  winds  force  the  plasma  to  move  up  and  down  the  magnetic  field  lines 
and  result  in  magnetic  field  induced  vertical  drifts.  The  general  trend  of  the  low  lat¬ 
itude  ( 5°N-5°S )  HWM  zonal  winds  (binned  every  hour  by  local  time  and  averaged 
over  all  Kp  and  F10.7)  is  westward  during  the  day  and  eastward  at  night  with  only 
slight  variations  in  magnitude  and  time  of  wind  direction  switch  (Figure  23).  Liu 
et  al.  (2006),  also  evaluated  the  HWM  zonal  winds  for  different  levels  of  geomagnetic 
activity  and  F10.7  values,  however,  the  variation  in  the  results  remains  small.  The 
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Ne  Percent  Increase  at  0°E  with  a  Doubled  0+/0  Collision  Frequency  for  Solar  Minimum  December  Solstice 
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increase  the  upper  altitude  plasma 


HWM  is  a  climatological  model  and  cannot  accurately  represent  the  global  thermo¬ 
spheric  wind  pattern  on  a  day  to  day  basis  and  is  instead  an  averaged  representation 
leading  to  uncertainties  in  the  model  output. 


Figure  23.  HWM  derived  zonal  thermospheric  winds  averaged  over  Kp  and  I'\ 0.7  at 
400fcm  for  all  three  geophysical  conditions.  The  winds  are  westward  (negative  val¬ 
ues)  during  the  day  and  eastward  (positive  values)  at  night  with  very  little  seasonal 
dependence  (Adapted  from  Liu  et  al.  (2006)) 


In  its  simplest  version,  the  earth’s  magnetic  field  can  be  represented  by  a  tilted 
and  offset  dipole  with  respect  to  earth’s  axis  of  rotation.  The  offset  dipole  field  results 
in  a  non-perpendicular  intersection  of  a  plane  of  magnetic  field  fines  and  the  earth’s 
surface  over  certain  geographic  locations.  Figure  24  shows  a  simple  schematic  of 
this  effect.  The  offset  dipole  field  and  the  daytime  westward  thermospheric  winds 
cause  upward  plasma  drifts  near  45°  E  (Madagascar)  and  downward  plasma  drifts 
near  225 °E  (Southeast  Pacific)  during  the  daytime;  these  drifts  are  reversed  at  night. 

In  addition  to  the  offset  dipole,  the  east-west  tilt  of  the  field  lines  also  leads 
to  magnetic  field  induced  vertical  drifts.  Figure  25  shows  the  HWM-derived  zonal 
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Figure  24.  Schematic  of  an  offset  dipole  magnetic  field  and  the  daytime  low-latitude 
thermospheric  winds  as  viewed  from  above  the  northern  geographic  pole.  The  zonal 
winds  force  the  ions  along  the  field  lines  resulting  in  upward  drifts  near  Madagascar 
and  downward  drifts  in  the  Southeast  Pacific 


winds  for  the  three  geophysical  conditions  used  in  this  project.  The  times  chosen 
(1000UT  and  1800UT  during  solar  minimum  and  solar  medium,  and  0000UT  and 
0500UT  during  solar  maximum)  correspond  to  the  time  of  the  maximum  change  in 
TEC  and  Ne.  Eastward  winds  are  depicted  as  positive  values  and  westward  winds  are 
negative  values;  in  general,  the  daytime  winds  are  westward  and  the  nighttime  winds 
are  eastward  as  was  seen  in  Figure  23.  The  black  lines  are  magnetic  field  lines  (and 
the  geomagnetic  equator)  as  calculated  by  the  IGRF.  The  IGRF  is  a  numerical  model 
that  calculates  earth’s  magnetic  field  based  on  observations  collected  by  satellites  at 
magnetic  observatories  and  during  magnetic  surveys  ( Finlay  et  al. ,  2010).  The  lowest 
order  solution  of  the  IGRF  is  a  dipole  field.  Because  of  this  magnetic  field  structure, 
the  zonal  wind  has  components  both  parallel  and  perpendicular  to  the  field  lines.  The 
parallel  component  will  cause  the  ions  to  move  along  the  field  lines;  if  the  parallel 
component  of  the  wind  is  toward  the  equator,  upward  drifts  will  occur.  Downward 
drifts  will  occur  with  a  parallel  component  that  is  poleward. 
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For  example,  in  the  upper  left  plot  in  Figure  25,  the  westward  flow  over  South 
Africa  causes  the  ions  to  move  along  the  field  line  toward  the  equator,  which  results  in 
an  upward  vertical  drift.  In  the  upper  right  plot,  the  westward  flow  over  the  Southeast 
Pacific  causes  the  ions  to  move  along  the  field  line  toward  the  pole  resulting  in  a 
downward  vertical  drift.  These  magnetic  field  induced  drifts  are  especially  important 
for  the  solar  maximum  June  solstice  case  (bottom  two  plots)  in  the  region  of  the 
South  Atlantic  Anomaly  (SAA)  which  will  be  discussed  further  below. 

For  this  comparison,  the  IPM  was  first  run  with  the  HWM-derived  zonal  winds 
and  then  with  the  zonal  winds  set  to  zero  to  evaluate  the  effects  of  an  uncertain  zonal 
wind  field  on  the  magnetic  field  induced  vertical  drifts.  Although  all  of  the  magnetic 
field  induced  drifts  are  shut  off  when  the  HWM  zonal  winds  are  forced  to  zero  in 
the  comparison  runs,  the  E  x  B  drifts  are  unaffected.  In  the  IPM,  each  parameter 
is  determined  using  separate  climatological  models  and  the  zonal  wind  model  is  not 
coupled  to  the  electric  field  model.  Therefore,  forcing  the  zonal  winds  to  zero  does 
not  affect  the  E  x  B  drift  that  is  calculated  within  the  model. 

The  most  significant  results  from  this  comparison  occur  during  the  day  near  Mada¬ 
gascar  (45° E)  and  in  the  Southeast  Pacific  Ocean  near  270° E.  The  largest  changes  in 
TEC  occur  over  Madagascar  at  1000UT  (1300LT)  and  over  the  Southeast  Pacific  at 
1800UT  (1100LT)  during  both  solar  minimum  December  solstice  and  solar  medium 
fall  equinox  conditions.  The  largest  changes  in  TEC  during  solar  maximum  June 
solstice  occur  in  the  Eastern  Pacific  and  the  South  Atlantic.  All  of  the  changes  are 
directly  related  to  the  complicated  structure  of  the  IGRF  derived  magnetic  field  and 
to  the  changes  in  the  magnetic  field  induced  vertical  drifts  of  the  plasma. 

When  the  zonal  winds  are  set  to  zero,  the  TEC  near  Madagascar  decreases  a 
maximum  of  45%  (12  TECU)  during  solar  minimum  and  50%  (40  TECU)  during 
solar  medium  (Figures  26  and  27).  This  occurs  because  the  magnetic  field  induced 
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300 km  HWM  Zonal  Winds  for  All  Three  Geophysical  Cases 
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Figure  25.  HWM  zonal  winds  at  300 km  for  all  three  geophysical  cases  at  1000UT  (left) 
and  1800UT  (right)  for  solar  minimum  (top)  and  solar  medium  (middle)  and  at  0000UT 
(left)  and  0500UT  (right)  for  solar  maximum  (bottom).  Positive  values  are  eastward 
and  negative  values  are  westward.  The  black  lines  are  magnetic  field  lines  as  calculated 
by  the  IGRF.  The  combination  of  the  zonal  winds  and  the  east-west  tilt  of  the  field 
lines  creates  vertical  plasma  drifts 
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upward  plasma  drift  is  stopped  during  the  day  keeping  the  plasma  at  lower  altitudes 
where  recombination  is  faster.  The  electron  density  decrease  at  45°  E  is  also  shown 
in  the  latitude  versus  altitude  plots  in  Figures  26  and  27.  Densities  decrease  above 
the  F2  peak  due  to  the  lack  of  the  upward  drift  and  increase  slightly  below  the  peak 
as  more  plasma  is  kept  at  lower  altitudes. 

TEC  over  the  Southeast  Pacific  increases  up  to  85%  (10  TECU)  during  solar 
minimum  and  57%  (16  TECU)  during  solar  medium  as  seen  in  Figures  28  and  29.  In 
this  region,  setting  the  zonal  winds  to  zero  shuts  off  the  downward  plasma  drift  keeping 
the  plasma  at  higher  altitudes  where  recombination  rates  are  slower.  Once  again,  the 
latitude  versus  altitude  plots  show  additional  information  about  the  altitude  at  which 
the  changes  in  electron  density  are  occurring.  The  bottom  right  difference  plot  shows 
an  increase  in  electron  densities  above  the  F2  peak  and  a  decrease  in  electron  densities 
below  the  peak.  This  is  the  opposite  from  what  is  occurring  over  Madagascar  because 
in  this  region  the  magnetic  field  induced  daytime  drift  is  opposite  that  of  Madagascar. 

The  most  significant  changes  in  TEC  that  occurred  during  solar  maximum  (Fig¬ 
ure  30)  were  different  than  for  the  other  two  cases.  Although  the  decrease  near  Mada¬ 
gascar  occurs  for  this  case  with  a  maximum  decrease  of  50%  (29  TECU)  at  0600UT, 
the  increase  in  the  Southeast  Pacific  does  not  occur.  The  reason  is  partly  due  to 
encroaching  eastward  winds  in  the  HWM  at  the  time  when  the  maximum  changes 
are  expected.  This  encroaching  eastward  wind  occurs  due  to  the  short  amount  of 
daylight  during  this  period  (southern  hemisphere  winter) .  The  TEC  over  the  Eastern 
Pacific  decreases  a  maximum  of  50  TECU  (90%)  and  the  TEC  over  the  South  Atlantic 
increases  75  TECU  (400%)  (Figure  30).  Both  of  these  changes  occur  at  night  (around 
2200LT)  when  the  eastward  winds  are  strongest  over  these  locations  (Figure  25). 

The  decrease  over  the  Eastern  Pacific  is  a  result  of  the  interaction  of  the  magnetic 
field  and  the  eastward  zonal  winds  creating  an  upward  drift  in  this  region.  Shutting 
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off  the  winds  stops  the  upward  drift  resulting  in  decreased  plasma  densities.  The 
increase  over  the  South  Atlantic  is  a  result  of  a  complicated  magnetic  held.  The 
magnetic  held  vertical  drifts  create  regions  of  converging  and  diverging  winds  which 
results  in  enhancements  and  depletions  in  the  plasma  density.  For  this  case,  the  result 
of  setting  the  zonal  winds  to  zero  is  a  signihcant  increase  in  the  equatorial  anomalies 
in  the  South  Atlantic  (upper  right  hand  plot,  Figure  30). 

4.5  Daytime  Production 

Daytime  production  mechanisms  include  photoionization,  chemical  reactions,  sec¬ 
ondary  electron  production,  and  impact  ionization  as  discussed  in  Chapter  IF  Sec¬ 
ondary  electron  production  is  an  additional  important  production  mechanism  that  is 
difficult  to  accurately  model.  This  leads  to  uncertainties  in  the  magnitude  of  the  total 
daytime  production  in  the  ionosphere.  Previous  studies  found  that  daytime  produc¬ 
tion  in  a  physics-based  model  needed  to  be  increased  by  a  factor  of  1.8  (the  ‘primo 
fix’)  to  account  for  the  contribution  of  secondary  electron  production.  More  recent 
studies  by  Zhu  et  al.  (2006)  using  a  larger  set  of  observational  data  have  suggested 
that  this  value  is  too  high  and  that  realistic  production  levels  are  obtained  when 
the  production  multiplication  factor  varies  linearly  as  a  function  of  F10.7  as  shown  in 
Figure  31. 

For  this  comparison,  the  baseline  runs  of  the  IPM  were  done  with  the  original  1.8 
multiplication  factor  while  the  comparison  runs  were  done  with  a  multiplication  factor 
that  varied  linearly  as  a  function  of  F10.7  (Figure  31).  The  maximum  change  to  the 
production  occurred  during  solar  maximum  (FW7  =  220)  with  the  production  multi¬ 
plication  factor  reduced  from  1.8  to  1.2 — a  decrease  of  33%  in  overall  production.  The 
minimum  change  occurred  during  solar  minimum  (F10.7  =  80)  with  a  new  production 
multiplication  factor  of  1.5,  which  is  a  decrease  of  17%  in  overall  production. 
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Figure  31.  The  decrease  in  daytime  production  multiplication  factor  varies  linearly  as 
a  function  of  F10.7.  The  largest  change  to  the  production  factor  occurred  for  high  F10.7 
values  and  the  smallest  change  occurred  for  low  F10.7  values 

Higher  solar  flux  during  solar  maximum  results  in  higher  production  due  to  both 
photoionization  as  well  as  to  secondary  electrons.  This  is  why  the  solar  maximum 
photoionization  production  is  only  increased  by  20%  (production  multiplication  factor 
of  1.2)  to  account  for  secondary  electron  production  while  during  solar  minimum  it 
is  increased  by  50%  (production  multiplication  factor  of  1.5).  For  example,  TEC 
during  solar  minimum  can  be  as  low  as  25  TECU;  increasing  this  by  50%  results 
in  a  12.5  TECU  increase.  Solar  maximum  TEC  is  as  large  as  100  TECU;  the  20% 
solar  maximum  increase  results  in  a  20  TECU  increase  which  is  7.5  TECU  higher 
than  the  solar  minimum  increase.  Therefore  the  higher  number  of  energetic  electrons 
during  solar  maximum  are  still  being  accounted  for  with  the  linear  variation  of  the 
production  multiplication  factor. 

The  results  of  this  comparison  show  that  the  magnitude  of  the  maximum  decrease 
in  TEC  during  the  day  is  directly  proportional  to  the  change  in  the  production  mul¬ 
tiplication  factor.  The  smallest  change  in  TEC  occurred  during  solar  minimum  while 
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solar  maximum  had  the  largest  change  in  TEC.  Figure  32  shows  the  TEC  plots  for 
solar  minimum  December  solstice.  The  difference  plot  shows  a  7  TECU  maximum 
decrease  occurring  at  1400LT  in  the  equatorial  anomalies.  Figures  33  and  34  show 
the  same  plots  for  both  solar  medium  fall  equinox  and  solar  maximum  June  solstice, 
respectively.  The  maximum  decrease  during  solar  medium  fall  equinox  was  30  TECU 
while  during  solar  maximum  June  solstice  the  maximum  decrease  was  50  TECU.  Ad¬ 
ditionally,  all  three  geophysical  cases  had  very  little  or  no  change  in  TEC  at  night 
when  production  due  to  photoionization  and  secondary  electron  production  decrease. 
Table  5  gives  a  summary  of  the  maximum  change  in  TEC  for  all  three  geophysical 
cases. 


Table  5.  TEC  decrease  with  decreased  daytime  production.  The  maximum  change  in 
TEC  was  directly  proportional  to  the  change  in  the  production  multiplication  factor 


Solar  Minimum 

Solar  Medium 

Solar  Maximum 

Original  Production  Factor 

1.8 

1.8 

1.8 

New  Production  Factor 

1.5 

1.35 

1.2 

Production  Factor  Decrease 

17% 

25% 

33% 

TECU  Decrease 

7 

30 

50 

TEC  Percent  Decrease 

17% 

27% 

33% 

4.6  Nighttime  E  x  B  Drift 

Photoionization  decreases  significantly  at  night  and  results  in  a  significant  decay  of 
the  ionosphere.  During  this  time,  the  ionosphere  is  primarily  controlled  by  loss  mech¬ 
anisms  and  plasma  transport.  At  low  latitudes,  the  nighttime  ExB  drift  is  downward 
due  to  the  westward  electric  field.  This  drift  brings  the  ions  to  lower  altitudes  where 
recombination  rates  are  higher  resulting  in  the  depletion  of  the  ionosphere.  Down¬ 
ward  drifts  that  are  too  fast  will  result  in  nighttime  densities  that  are  too  low  while 
drifts  that  are  too  slow  will  result  in  densities  that  are  too  high  due  to  the  higher 
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TEC  with  Decreased  Daytime  Production  for  Solar  Minimum  December  Solstice 
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Figure  32.  TEC  with  decreased  daytime  production  for  case  1  (solar  minimum  Decem¬ 
ber  solstice)  at  OOOOUT.  The  top  left  plot  is  with  a  1.8  multiplication  factor  for  daytime 
production,  the  top  right  plot  is  with  a  1.5  multiplication  factor,  and  the  bottom  plot 
is  the  difference.  The  17%  decrease  in  the  production  multiplication  factor  resulted  in 
a  maximum  decrease  of  7  TECU  (17%)  in  the  equatorial  anomalies 
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TEC  with  Decreased  Daytime  Production  for  Solar  Medium  Fall  Equinox 


Production  Coefficient  =  1.8 
Solar  Med  Fall  Equinox  (1998  267)  OOOOUT 


Production  Coefficient  =  1 .35 


Figure  33.  TEC  with  decreased  daytime  production  for  case  2  (solar  medium  fall 
equinox)  at  OOOOUT.  The  top  left  plot  is  with  a  1.8  multiplication  factor  for  daytime 
production,  the  top  right  plot  is  with  a  1.35  multiplication  factor,  and  the  bottom  plot 
is  the  difference.  The  25%  decrease  in  the  production  multiplication  factor  resulted  in 
a  maximum  decrease  of  30  TECU  (27%)  in  the  equatorial  anomalies 
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TEC  with  Decreased  Daytime  Production  for  Solar  Maximum  June  Solstice 


Production  Coefficient  =  1.8 
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Figure  34.  TEC  with  decreased  daytime  production  for  case  3  (solar  maximum  June 
solstice)  at  0000UT.  The  top  left  plot  is  with  a  1.8  multiplication  factor  for  daytime 
production,  the  top  right  plot  is  with  a  1.2  multiplication  factor,  and  the  bottom  plot 
is  the  difference.  The  33%  decrease  in  the  production  multiplication  factor  resulted  in 
a  maximum  decrease  of  50  TECU  (33%)  in  the  equatorial  anomalies 
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recombination  at  low  altitudes.  Slowing  the  downward  nighttime  E  x  B  drift  allows 
the  ions  to  remain  at  higher  altitudes  where  recombination  rates  are  slower.  In  this 
comparison,  the  downward  nighttime  E  x  B  drift  is  decreased  linearly  as  a  function 
of  F10.7  (Figure  35)  to  investigate  the  resulting  changes  to  the  nighttime  densities. 


Downward  E  x  B  Drift  Decrease 


Figure  35.  Decrease  in  downward  E  x  B  drift  varies  linearly  as  a  function  of  F10.7.  The 
largest  change  occurred  for  low  and  the  smallest  change  occurred  for  large  F10.7  values 


For  all  three  geophysical  cases,  the  largest  absolute  increase  in  nighttime  TEC 
occurs  around  0100LT-0200LT  while  the  largest  percent  increase  in  TEC  occurs  at 
0500LT  (due  to  the  much  smaller  densities  at  this  time) .  The  largest  changes  for  solar 
minimum  December  solstice  and  solar  medium  fall  equinox  occurred  in  the  region  of 
the  SAA  (Figures  36  and  37).  The  magnetic  field  is  weakest  in  this  region  resulting 
in  a  fast  initial  E  x  B  drift  (ue  oc  E/B).  Decreasing  this  fast  drift  results  in  a  larger 
overall  change  than  in  the  regions  where  the  drift  is  slower.  For  solar  maximum  June 
solstice,  the  largest  changes  occurred  over  and  just  to  the  west  of  South  America 
(Figure  38)  due  to  the  magnetic  field  configuration  and  neutral  wind  pattern  in  this 
region. 
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For  the  solar  minimum  December  solstice  case,  the  E  x  B  drift  was  decreased 
by  70%  resulting  in  a  maximum  8  TECU  increase  (480%  increase)  over  the  SAA  as 
seen  in  the  right  hand  plot  of  Figure  36.  A  slight  decrease  in  TEC  (less  than  1.5 
TECU)  also  occurs  during  the  day  due  to  the  residual  effects  propagating  into  the 
daytime  sector.  The  bottom  plots  in  Figure  36  show  the  45° E  altitude  structure  of 
the  electron  density  at  0500LT  (0300UT).  The  increase  in  low  latitude  densities  can 
be  seen  in  the  right  hand  plot.  This  plot  shows  a  maximum  700%  increase  in  electron 
density  occurring  over  the  geomagnetic  equator  at  approximately  300 Arm.  This  plot 
also  shows  that  the  changes  occur  only  in  the  low  latitudes — 30°  AV-200  S  geographic 
latitude  as  this  is  where  the  E  x  B  drift  is  most  significant.  Also  of  interest  is  the 
decrease  in  electron  densities  around  200Arm  in  the  southern  hemisphere.  The  decrease 
is  occurring  because  the  plasma  is  not  allowed  to  move  to  lower  altitudes  to  replenish 
the  depleted  plasma.  The  effect  is  stronger  in  the  southern  (summer)  hemisphere  due 
to  the  winter  anomaly. 

The  same  general  results  were  obtained  for  the  solar  medium  fall  equinox  case 
as  for  solar  minimum.  The  maximum  increase  in  TEC  occurred  over  the  SAA;  this 
increase  can  be  easily  seen  in  the  top  row  of  plots  in  Figure  37.  The  right  hand  plot 
shows  that  the  maximum  increase  in  TEC  was  almost  40  TECU,  corresponding  to 
a  630%  increase.  The  small  decrease  during  the  day  can  also  be  seen  in  this  plot. 
The  plots  showing  the  altitude  structure  of  the  electron  density  (Figure  37,  bottom 
row)  show  the  same  overall  changes  as  for  the  previous  case.  Increases  in  low  latitude 
electron  density  of  up  to  900%  result  from  decreasing  the  downward  nighttime  E  x  B 
drift  by  approximately  55%. 

The  results  for  the  solar  maximum  June  solstice  case  are  slightly  different  than 
the  previous  two  cases  due  to  the  significant  increase  in  photoionization  during  this 
period.  In  addition,  this  period  also  had  the  smallest  decrease  (30%)  to  the  downward 
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nighttime  E  x  B  which  results  in  smaller  changes  to  the  TEC  and  Ne.  The  top  right 
hand  plot  in  Figure  38  shows  that  the  maximum  increase  in  TEC  is  30  TECU  (160% 
increase)  while  the  bottom  right  hand  plot  shows  the  maximum  increase  in  electron 
densities  with  altitude  is  250%.  The  effect  of  the  winter  anomaly  can  also  be  seen  in 
this  plot  with  the  northern  (summer)  hemisphere  having  a  larger  decrease  in  electron 
densities  below  the  peak. 

Another  interesting  way  to  examine  the  results  of  this  comparison  is  to  look  at 
the  same  local  time  for  all  locations.  Figure  39  shows  both  the  absolute  difference 
(left)  and  percent  increase  (right)  in  TEC  at  0500LT  for  all  three  cases.  For  the  solar 
minimum  and  solar  medium  cases,  this  plot  shows  that  the  largest  changes  occur  over 
the  SAA  as  previously  stated.  However,  these  plots  also  show  that  the  location  of 
the  largest  changes  moves  from  the  SAA  during  solar  minimum  December  solstice  to 
the  East  Asia/Pacific  Ocean  region  during  solar  maximum  June  solstice.  This  change 
can  be  attributed  to  the  configuration  of  the  magnetic  field  and  the  differences  in 
field-aligned  drifts  that  result  from  the  changing  thermospheric  winds  with  season. 

4.7  Four- Wave  Tidal  Influence 

A  ‘longitudinal  wave  number  four  pattern’  of  enhanced  TEC  was  shown  to  exist  in 
the  low  latitude  equatorial  anomalies  by  Immel  et  al.  (2006).  This  wave  number  four 
pattern  is  created  by  longitudinal  variations  in  atmospheric  tides  in  the  troposphere 
that  propagate  upward  into  the  E  region  ionosphere  resulting  in  a  corresponding 
modulation  of  the  zonal  electric  field  and  the  vertical  plasma  drifts.  These  variations 
manifest  in  the  ionosphere  as  a  wave-like  pattern  of  TEC  enhancements  and  deple¬ 
tions  at  equatorial  and  low  latitudes.  Scherliess  et  al.  (2008)  studied  the  longitudinal 
variability  of  the  low  latitude  TEC  using  TOPEX  observations  and  showed  the  lo¬ 
cal  time,  seasonal,  solar  cycle,  and  geomagnetic  activity  dependence  of  this  pattern. 
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TEC  and  Ne  with  Decreased  Downward  E  x  B  Drift  for  Solar  Maximum  June  Solstice 
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TEC  with  Decreased  Downward  E  x  B  Drift  for  All  Three  Geophysical  Cases 
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Figure  39.  TEC  increase  with  decreased  downward  E  x  B  drift  for  all  three  geophysical 
cases  (solar  minimum  December  solstice  -  top,  solar  medium  fall  equinox  -  middle,  and 
solar  maximum  June  solstice  -  bottom).  The  plots  show  the  changes  in  TEC  at  0500LT 
for  all  geographic  locations.  The  area  of  the  maximum  change  shifts  from  the  SAA 
to  the  Pacific  Ocean  from  case  1  to  case  3.  Note  that  the  scales  are  different  for  each 
geophysical  case 
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It  was  found  that  the  wave  number  four  pattern  is  created  during  the  day  during 
equinox  and  June  solstice  and  is  present  but  not  well-defined  during  December  sol¬ 
stice  ( Scherliess  et  al,  2008).  It  was  also  found  that  the  pattern  is  not  dependent  on 
solar  cycle  or  geomagnetic  activity. 

Because  the  wave  number  four  pattern  is  created  by  tides  in  the  troposphere,  this 
forcing  can  only  be  simulated  in  an  ionospheric  model.  This  was  done  in  the  IPM  by 
modulating  the  vertical  E  x  B  drift  as  a  function  of  longitude  based  on  previous  work 
done  by  USU  ( Schunk ,  2010).  This  modulation  increases  the  E  x  B  drift  at  certain 
longitudes  and  decreases  it  at  other  longitudes. 

Including  tidal  forcing  in  this  manner  resulted  in  a  distinctive  pattern  of  en¬ 
hancements  and  depletions  of  TEC  and  400 km  Ne  in  the  equatorial  anomalies  (Fig¬ 
ures  40  and  41).  The  enhancements  are  centered  at  longitudes  of  15E,  110E,  200E, 
and  290E  and  the  depletions  centered  at  longitudes  of  60E,  155E,  245E,  and  335E 
this  pattern  does  not  migrate  with  time.  Each  enhancement  and  depletion  spans 
approximately  40°  of  longitude  and  is  symmetric  in  location  (though  not  magnitude) 
about  the  geometric  equator  for  all  three  geophysical  cases. 

Figures  40  and  41  show  the  TEC  and  400 km  Ne,  respectively,  for  solar  minimum 
December  solstice.  Both  figures  are  at  1300LT,  which  is  when  the  largest  changes 
in  each  parameter  occur.  The  top  left  plot  in  each  figure  shows  the  model  output 
without  the  E  x  B  drift  modulation  while  the  top  right  plot  is  the  model  output 
with  the  tidal  forcing  included.  Although  these  two  plots  look  nearly  identical,  sub¬ 
tracting  one  grid  from  the  other  shows  the  enhancements  and  depletions  that  result 
from  the  modulation  (bottom  plot).  For  this  case,  the  TEC  increases  a  maximum 
of  23%  over  Northern  Africa  and  decreases  a  maximum  of  14%  over  the  Northeast 
Pacific.  Changes  in  the  400A:m  Ne  are  similar  with  a  maximum  increase  of  44%  and 
a  decrease  of  24%  in  the  same  locations.  Of  interest  is  that  the  northern  hemisphere 
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enhancements  and  depletions  are  slightly  greater  than  the  southern  hemisphere  ones. 
This  is  again  due  to  the  winter  anomaly  and  the  resulting  asymmetries  in  electron 
densities. 

Because  the  differences  in  the  plots  of  TEC  and  400 km  Ne  without  tidal  forcing 
and  with  tidal  forcing  (e.g.  between  Figures  40  and  41)  are  difficult  to  distinguish, 
these  plots  are  not  shown  for  the  solar  medium  and  solar  maximum  cases.  Instead, 
Figure  42  shows  only  the  plots  of  the  percent  change  in  TEC  and  400 km  Ne  for  both 
cases.  As  mentioned  previously  and  seen  clearly  in  these  plots,  the  pattern  of  en¬ 
hancements  and  depletions  is  similar  across  solar  cycle  and  season  with  only  a  change 
in  the  magnitude  of  the  enhancements  and  depletions.  Table  6  lists  the  maximum 
changes  in  TEC  and  400 km  Ne  for  all  three  cases.  The  smallest  percent  change  occurs 
during  solar  maximum  June  solstice  (bottom  row  of  Figure  42),  which  is  due  to  the 
higher  overall  plasma  densities  present  during  these  conditions.  Again,  the  effects 
of  the  winter  anomaly  can  be  seen  in  the  asymmetric  magnitudes  in  the  June  sol¬ 
stice  case  (bottom  row)  with  the  southern  hemisphere  having  larger  changes  than  the 
northern  hemisphere;  for  the  fall  equinox  case,  the  magnitudes  of  the  enhancements 
and  depletions  are  fairly  symmetric  across  the  geomagnetic  equator. 


Table  6.  Maximum  percent  change  in  TEC  and  400 km  X,  at  1300LT.  The  largest 
changes  occur  during  solar  minimum  December  solstice  while  the  smallest  changes 
occur  during  solar  maximum  June  solstice 


TEC 

400km  Ne 

Maximum 
%  Increase 

Maximum 
%  Decrease 

Maximum 
%  Increase 

Maximum 
%  Decrease 

Case  1 

23 

14 

44 

24 

Case  2 

21 

16 

34 

28 

Case  3 

13 

11 

19 

18 
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TEC  with  Tidal  Forcing  Included  for  Solar  Minimum  December  Solstice 


Without  Tidal  Forcing 

Solar  Min  Dec  Solstice  (1995  358)  1300LT 


TEC  Units 


With  Tidal  Forcing 

Solar  Min  Dec  Solstice  (1995  358)  1300LT 


TEC  Units 


TEC  Percent  Increase 
Solar  Min  Dec  Solstice  (1995  358)  1300LT 


TEC  Percent  Increase 


Figure  40.  TEC  with  tidal  forcing  for  case  1  (solar  minimum  December  solstice).  The 
plots  are  for  a  fixed  local  time  (1300LT).  The  top  left  plot  is  the  default  run,  the 
top  right  plot  is  with  tidal  forcing  included,  and  the  bottom  plot  is  the  difference. 
Modulating  the  E  x  B  drift  results  in  enhancements  and  depletions  in  the  low  latitude 
TEC 
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Ne  at  400 km  with  Tidal  Forcing  Included  for  Solar  Minimum  December  Solstice 
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Figure  41.  Ne  at  400 km  with  tidal  forcing  for  case  1  (solar  minimum  December  solstice). 
The  plots  are  for  a  fixed  local  time  (1300LT).  The  top  left  plot  is  the  default  run,  the 
top  right  plot  is  with  tidal  forcing  included,  and  the  bottom  plot  is  the  difference. 
Modulating  the  E  x  B  drift  results  in  enhancements  and  depletions  in  the  low  latitude 

400 km  Ne 
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TEC  and  400 km  Ne  Percent  Increase  with  Tidal  Forcing  Included 
for  Two  Geophysical  Cases 
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Solar  Max  June  Solstice  (2001  172)  1300LT 
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Figure  42.  TEC  and  400 km  Ne  with  tidal  forcing  for  case  2  and  3.  The  left  hand  plots 
are  the  percent  increase  in  TEC  and  the  right  hand  plots  are  the  percent  increase 
in  400 km  Ne  for  solar  medium  fall  equinox  (top)  and  solar  maximum  June  solstice 
(bottom).  The  plots  are  for  a  fixed  local  time  (1300LT).  Modulating  the  E  x  B  drift 
results  in  enhancements  and  depletions  in  the  low  latitude  TEC  and  400 km  Ne.  Note 
that  the  scales  are  different  for  each  geophysical  case 


74 


V.  Conclusions  and  Recommendations 


5.1  Chapter  Overview 

This  chapter  is  broken  into  two  sections.  The  first  section  summarizes  the  analysis 
of  the  physical  parameters  in  the  IPM  and  reiterates  the  most  significant  changes  that 
occurred  due  to  changing  these  parameters.  Recommendations  for  future  research  are 
given  in  the  second  section  of  this  chapter. 

5.2  Conclusions 

Understanding  the  effects  of  adjusting  parameters  in  a  physics-based  model  is 
critical  in  ensuring  the  accuracy  of  the  model.  This  study  analyzed  model  performance 
by  varying  the  0+/0  collision  frequency,  downward  nighttime  E  x  B  drifts,  daytime 
production,  zonal  winds,  and  tidal  forcing.  It  was  found  that  when  the  physical 
parameters  were  adjusted  in  the  IPM,  the  resulting  changes  were  nonlinear  and  varied 
with  location,  time,  and  solar  conditions.  The  changes  ranged  from  very  slight  changes 
in  the  model  output  to  a  630%  increase  in  TEC  and  a  900%  increase  in  Ne.  Each 
parameter  had  its  own  unique  effect  on  the  model  output;  the  most  significant  of 
these  changes  are  summarized  below. 

The  most  significant  effect  of  doubling  the  0+/0  collision  frequency  was  an  in¬ 
crease  in  the  mid  latitude  nighttime  peak  electron  density  ranging  from  30-180%. 
The  changes  were  not  consistent  across  solar  cycle,  season,  or  location  and  so  no 
generalities  could  be  made  as  to  when  the  smallest  and  largest  changes  occured.  Also 
significant  were  decreases  in  electron  densities  below  the  peak  with  a  maximum  80% 
decrease.  This  decrease  occured  because  the  doubled  0+/0  collision  frequency  drove 
the  plasma  to  higher  altitudes  and  the  plasma  depleted  by  recombination  at  low 
altitudes  was  not  replenished. 
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It  was  found  that  setting  the  zonal  winds  to  zero  resulted  in  an  increase  in  TEC 
of  up  to  400%  over  the  Southeast  Pacific  Ocean  with  the  largest  change  occurring 
during  solar  maximum  June  solstice  at  1100LT.  Near  Madagascar  on  the  opposite 
side  of  the  globe,  the  TEC  decreased  with  a  maximum  change  of  50%  during  solar 
medium  fall  equinox  at  1300LT.  The  changes  in  TEC  were  caused  by  changes  to  the 
magnetic  field  induced  vertical  drifts  created  by  the  combination  of  zonal  winds  and 
the  earth’s  magnetic  field. 

The  most  straightforward  results  occurred  when  the  production  multiplication  fac¬ 
tor  was  scaled  to  correctly  account  for  secondary  electron  production.  The  changes  to 
Ne  and  TEC  were  found  to  be  directly  proportional  to  the  change  in  the  production 
factor.  The  solar  maximum  production  multiplication  factor  was  decreased  by  33% 
with  a  resulting  33%  decrease  in  low  latitude  TEC.  A  production  multiplication  fac¬ 
tor  decrease  of  25%  and  17%  during  solar  medium  and  minimum,  respectively,  also 
resulted  in  directly  proportional  changes  to  the  TEC  in  the  equatorial  anomalies. 

The  most  significant  results  of  decreasing  the  nighttime  downward  E  x  B  drift 
occurred  during  solar  medium  fall  equinox.  TEC  values  over  the  SAA  increased  630% 
while  300 km  Ne  increased  900%.  The  solar  minimum  December  solstice  case  had  only 
slightly  smaller  changes,  while  the  solar  maximum  June  solstice  case  had  increases 
of  160%  and  250%  in  TEC  and  300 km  Ne,  respectively.  These  large  changes  in  the 
model  output  show  the  importance  of  correctly  specifying  the  physical  parameters. 

For  the  final  adjustment,  tidal  forcing  was  included  in  the  model  by  modulating 
the  E  x  B  drift  and  was  found  to  reproduce  the  four-wave  pattern  of  enhanced 
TEC  at  low  latitudes.  Low  latitude  TEC  and  400 km  Ne  increased  13-23%  and  19- 
44%,  respectively,  at  longitudes  centered  at  15 °E,  110° A,  200° E,  and  290° E  while 
decreasing  11-16%  (TEC)  and  18-28%  (400/cm  Ne)  at  longitudes  centered  at  60° E, 
155 °E,  245 °E,  and  335° E.  Because  the  same  modification  was  applied  to  all  three 
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cases,  the  largest  relative  changes  occured  during  solar  minimum  and  the  smallest 
relative  changes  occurred  during  solar  maximum. 

Except  for  the  adjustment  to  the  daytime  production  multiplication  factor,  none 
of  the  adjustments  produced  linear  or  straightforward  changes  in  the  model  output. 
Not  only  did  the  changes  vary  as  a  function  of  time  and  solar  condition,  but  they 
also  varied  with  location.  These  nonlinear  variations  with  location  resulted  from 
the  earth’s  complex  magnetic  field  configuration  and  the  zonal  wind  distribution  in 
addition  to  many  other  factors.  The  nonlinearity  of  these  relationships  makes  it 
difficult  to  identify  which  parameters  may  be  causing  erroneous  model  output. 

5.3  Recommendations  for  Future  Work 

Two  avenues  exist  for  possible  future  work  stemming  from  the  results  of  this 
project.  The  first  is  to  improve  the  operational  capability  of  the  IPM  by  comparing 
the  model  output  to  known  data  and  adjusting  the  parameters  until  the  results  are 
better  correlated.  The  second  avenue  is  to  investigate  some  of  the  physical  parameters 
more  in-depth  for  a  larger  set  of  geophysical  conditions. 

While  this  project  provided  useful  information  on  the  range  of  output  that  oc¬ 
curred  when  certain  physical  parameters  were  adjusted  in  the  model,  no  information 
was  given  as  to  the  accuracy  of  the  model  output,  either  adjusted  or  unadjusted.  In 
order  for  the  model  to  be  useful  operationally,  it  must  be  validated.  This  can  be  done 
by  comparing  the  model  output  either  to  another  previously  validated  physical  model 
(e.g.  the  IFM)  or  to  an  actual  data  set  (e.g.  the  TOPEX  database). 

Comparing  the  IPM  to  the  IFM  would  be  the  most  simple  and  straightforward 
way  to  perform  an  initial  validation  of  the  IPM.  The  physical  parameters  have  already 
been  optimized  in  the  IFM  and  using  these  optimized  values  for  the  parameters  in 
the  IPM  would  provide  an  ideal  starting  point  for  the  IPM  validation.  In  addition, 
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comparing  the  IPM  to  the  IFM  would  allow  a  wide  range  of  controlled  geophysical 
cases  to  be  examined.  Ideally,  this  comparison  would  serve  to  optimize  the  physical 
parameters  in  the  IPM  before  it  is  compared  to  actual  observational  data. 

Even  if  the  IPM  is  compared  to  another  physical  model,  the  IPM  still  needs  to 
be  compared  to  actual  data  to  ensure  the  accuracy  of  model.  Because  the  TOPEX 
TEC  database  covers  a  large  range  of  geophysical  conditions,  it  would  be  a  good  ob¬ 
servational  database  to  use  for  the  validation.  However,  there  are  many  other  sources 
of  data  that  can  be  used  to  validate  the  IPM  including  ionosondes,  GPS  TEC,  and 
the  Constellation  Observing  System  for  Meteorology,  Ionosphere,  and  Climate  (COS¬ 
MIC)  occultations.  If  an  IPM/IFM  type  comparison  is  not  done,  this  comparison 
would  serve  as  an  opportunity  to  optimize  the  physical  parameters  in  the  IPM. 

While  the  operational  capability  of  the  IPM  is  important,  this  research  also  has 
significant  implications  from  a  theoretical  viewpoint.  Understanding  in  more  depth 
how  adjusting  each  of  these  parameters,  specifically  the  0+/0  collision  frequency  and 
the  zonal  winds,  modifies  electron  densities  and  TEC  would  be  useful  in  gaining  a 
better  understanding  of  ionospheric  processes.  These  two  parameters  are  highlighted 
because  they  showed  the  greatest  nonlinearity  in  the  model  output  changes.  While 
this  project  provided  a  broad  overview  of  the  range  of  model  output  that  results  from 
the  uncertainty  in  these  physical  parameters,  these  two  parameters  could  be  looked 
at  in  more  detail.  For  example,  this  project  varied  the  level  of  solar  activity  and 
season  together;  different  relationships  may  be  present  if  only  the  seasons  are  varied 
and  the  level  of  solar  activity  is  kept  constant.  Additionally,  the  temporal  and  spatial 
resolution  could  be  increased  over  regions  where  the  model  output  is  changing  the 
most.  Focusing  the  research  on  just  these  two  parameters  for  more  geophysical  cases 
with  an  increased  temporal  and  spatial  resolution  may  show  additional  trends  and 
relations  in  the  nonlinearity  of  the  results  presented  in  this  paper. 
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Appendix  A.  Ionosphere  Forecast  Model 


The  Ionosphere  Forecast  Model  (IFM)  is  a  physics-based  model  that  calculates 
the  three-dimensional  time  dependent  electron  and  ion  density  distributions  from  90- 
1600 km  ( S chunk  et  al,  1997).  Density  distributions  are  given  for  four  major  ions 
(. NO+ ,  O2,  and  0+ )  in  the  E  region,  two  major  (0+  and  NO+ )  and  two  minor 

(N£  and  O2  )  ions  in  the  F  region,  and  the  ion  and  electron  temperatures  for  the  entire 
profile  [S chunk  et  al,  2004).  In  addition,  the  IFM  calculates  the  H+  densities  for 
the  F  region  and  the  topside  ionosphere.  A  numerical  solution  to  the  ion  continuity, 
momentum,  and  energy  equations  is  used  to  determine  the  ion  distributions.  Based 
on  the  quasi-neutral  state  of  the  ionosphere,  the  electron  density  profile  is  derived 
from  the  ion  density  information. 

In  the  E  region,  chemical  equilibrium  (Equation  2)  is  assumed  for  all  four  major 
ions  with  each  ion  having  comparable  density.  The  four  coupled  nonlinear  equations 
are  solved  using  a  Taylor  series  expansion  and  iterating  to  a  final  solution  (S chunk 
and  Sojka ,  1994).  In  the  F  region,  transport  is  important  and  the  two  major  ion 
densities  are  solved  using  numerical  techniques  for  the  coupled  continuity  and  mo¬ 
mentum  equations.  The  two  minor  ions  in  the  F  region  are  assumed  to  be  in  chemical 
equilibrium  and  their  densities  are  solved  using  Equation  2.  The  H+  densities  in  the 
topside  ionosphere  are  computed  by  comparing  the  concentrations  of  H+  and  0+  and 
by  using  a  combination  of  chemical  and  diffusive  equilibrium  relations. 

The  IFM  can  be  run  on  a  global,  regional,  or  local  scale  and  has  a  global  spatial 
resolution  of  3°  latitude  and  7.5°  longitude.  The  internal  vertical  resolution  of  the 
IFM  is  4 km,  but  the  output  resolution  is  variable  starting  at  4 km  in  the  E  region, 
increasing  to  16 km  in  the  F  region,  and  has  a  maximum  spacing  of  9 6km  at  1600£;m. 
The  magnetic  field  is  determined  using  a  simplified  version  of  the  International  Geo¬ 
magnetic  Reference  Field  (IGRF)  that  accounts  for  a  tilted  offset  dipole  field. 
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The  IFM  is  an  autonomous  model  that  uses  empirical  (statistical)  models  for  many 
of  the  required  inputs.  The  major  inputs  required  by  the  IFM  are  density  and  temper¬ 
ature  profiles  for  the  major  neutrals  (V2,  02,  and  O),  neutral  winds,  magnetospheric 
and  equatorial  electric  fields,  auroral  precipitation,  and  topside  heat  flow  ( Schunk 
and  Sojka ,  1994).  The  Mass  Spectrometer  and  Incoherent  Scatter  (MSIS)  model  is 
used  for  the  neutral  density  and  temperature  profiles.  This  empirical  model  provides 
these  parameters  based  on  a  given  set  of  geophysical  conditions  (solar  cycle,  season, 
geomagnetic  activity,  etc.).  The  zonal  and  meridional  neutral  winds  are  obtained 
from  the  Horizontal  Wind  Model  (HWM)  developed  by  Hedin  et  al.  (1991).  To  ac¬ 
count  for  the  magnetospheric  electric  field,  the  IFM  uses  the  Heppner  and  Maynard 
(1987)  model,  which  gives  electric  fields  in  the  polar  regions  as  a  function  of  Kp.  The 
dynamo  equatorial  electric  field  model  that  is  used  in  IFM  calculations  is  the  Fejer 
et  al.  (1999)  model.  Auroral  precipitation  accounted  for  in  the  IFM  using  the  model 
developed  by  Hardy  et  al  (1985).  Finally,  Schunk  et  al.  (1997)  discusses  how  the 
topside  heat  flow  is  modeled  in  the  IFM. 

Validation  of  the  IFM  resulted  in  minor  adjustments  to  certain  physical  parame¬ 
ters  in  the  model  due  to  discrepancies  between  model  output  and  observations.  These 
adjustments  included  modifications  to  the  0+/0  collision  frequency,  zonal  winds,  tidal 
forcing,  and  ion  production  rates.  The  0+/ O  collision  frequency  is  known  to  be  uncer¬ 
tain  up  to  a  factor  of  2.0;  correct  densities  were  obtained  in  the  IFM  with  a  coefficient 
of  1.0  ( Schunk ,  2010).  It  was  found  that  in  certain  geographic  regions,  most  notably 
around  Madagascar,  the  combination  of  the  magnetic  field  configuration  and  HWM 
derived  zonal  winds  led  to  large  upward  plasma  drifts  and  enhanced  electron  densities 
in  the  IFM  output.  Setting  the  zonal  winds  to  zero  corrected  these  density  enhance¬ 
ments.  The  third  adjustment  was  to  the  daytime  ion  production  rate.  This  rate 
was  decreased  as  a  linear  function  of  Fi0  7  ( Schunk ,  2010)  with  the  larger  decreases 
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occurring  during  solar  maximum  conditions.  This  adjustment  was  done  to  counter 
the  uncertainty  in  the  secondary  electron  production  calculation.  Another  adjust¬ 
ment  corrected  the  nighttime  equatorial  regions  with  TEC  values  that  were  too  low 
by  5-10  TEC  units.  USU  corrected  this  problem  by  decreasing  the  nighttime  E  x  B 
drift  and  tested  the  output  for  various  geophysical  conditions.  A  four- wave  pattern 
of  enhanced  TEC  is  known  to  exist  in  the  ionosphere  and  is  created  by  atmospheric 
tides  ( Scherliess  et  al. ,  2008).  The  IFM  cannot  reproduce  these  tides,  and  so  a  lon¬ 
gitudinal  modulation  of  the  E  x  B  drift  was  introduced  into  the  IFM  to  simulate 
the  tidal  forcing.  Adding  this  modulation  reproduced  the  known  pattern  of  TEC  en¬ 
hancements.  Finally,  the  high  latitude  portion  of  the  model  had  to  be  adjusted  due  to 
unrealistic  F  region  densities  near  South  America  for  specific  geophysical  conditions 
(winter,  solar  minimum,  06LT-12LT).  The  unique  combination  of  the  season,  winds, 
magnetic  field  configuration,  and  high-latitude  convection  in  this  region  formed  the 
enhanced  densities.  Table  7  summarizes  the  improvements  made  to  the  IFM. 


Table  7.  Improvements  to  the  Ionosphere  Forecast  Model 


Parameter 

Description 

0+/0  Collision  Frequency 

Scaled  the  0+/0  collision  frequency 

Zonal  Winds 

Set  zonal  winds  in  HWM  to  zero 

Tidal  Forcing 

Included  atmospheric  tidal  forcing  at  low  latitudes 

Daytime  Production 

Decreased  production  as  a  linear  function  of  Fjo.7 

Nighttime  Drifts 

— *  — * 

Decreased  nighttime  equatorial  E  x  B 

High  Latitudes 

Adjusted  electron  densities  near  South  America 

The  IFM  requires  as  input  the  day,  year,  UT  start  time,  number  of  loops  to  run, 
output  interval,  daily  and  90-day  average  F10.7  solar  flux,  three-hour  Kp  index,  daily 
Ap  index,  and  internal  and  external  resolution  settings.  If  a  previous  IFM  run  is  not 
available,  the  IFM  will  be  initiated  using  the  empirical  International  Reference  Iono¬ 
sphere  (IRI)  model.  This  is  not  ideal,  and  in  this  case  the  IFM  output  is  degraded  for 
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the  first  12  hours.  Therefore,  the  model  should  be  run  for  two  days — one  degraded 
day  and  one  good  day  of  output  ( Space  Environment  Corporation ,  2002).  Running 
the  model  for  longer  than  two  days  also  causes  invalid  output  as  some  of  the  adjust¬ 
ments  to  the  model  do  not  extend  past  day  two  ( Schunk ,  2010).  The  two-day  IFM 
run  takes  approximately  six  hours  and  the  output  is  in  a  scientific  binary  data  for¬ 
mat  as  a  Network  Common  Data  Form  (netCDF)  file.  The  IFM  output  includes  the 
one-dimensional  latitude,  longitude,  and  altitude  arrays;  the  two-dimensional  TEC, 
F2  and  E  region  peak  electron  density  and  height  arrays;  and  the  three-dimensional 
electron  density,  ion  density,  and  ion  and  electron  temperature  arrays.  More  infor¬ 
mation  on  the  options  available  in  the  IFM  and  how  to  run  the  IFM  is  available  in 
the  IFM  User  Manual  ( Space  Environment  Corporation,  2002). 
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Appendix  B.  Governing  Processes  in  the  Ionosphere 

2.1  Neutral  Atmosphere 

The  major  neutrals  in  the  ionosphere  are  molecular  nitrogen  (N2),  molecular  oxy¬ 
gen  (O2),  and  atomic  oxygen  (O),  with  minor  concentrations  of  argon  (Ar),  helium 
{He),  and  hydrogen  ( H ).  The  concentrations  of  each  of  these  species  varies  with 
height  and  can  be  approximated  using  a  hydrostatic  equilibrium  solution.  Hydro¬ 
static  equilibrium  occurs  when  the  forces  due  to  pressure  and  gravity  are  in  balance. 
Where  this  balance  occurs,  the  density  of  a  species  as  a  function  of  altitude  (z)  is 
given  as 

N8(z)  =  Ns(z0)  exp  (19) 

where  Ns  is  the  number  density,  z0  is  a  reference  altitude,  and  Hs  is  the  scale  height 
for  the  species  s  given  by 

=  khTs{z) 
msg(z ) 

where  Ay  is  Boltzmann’s  constant,  Ts  is  the  temperature,  ms  is  the  mass,  and  g  is 
gravity.  In  this  approximation,  the  concentration  of  the  species  falls  off  primarily 
as  a  function  of  mass.  Therefore,  the  concentrations  of  N2  and  O2  decrease  more 
rapidly  with  height  than  does  the  concentration  of  O.  Figure  43  shows  the  vertical 
variation  up  to  lOOOA’m  of  the  neutral  concentrations  during  solar  maximum  using 
this  assumption.  The  relative  abundance  of  each  neutral  compared  with  the  other 
neutrals  is  important  when  calculating  the  production  and  loss  of  ions. 

2.2  Production  and  Loss 

Photoionization  is  the  primary  production  mechanism  for  creating  ion-electron 
pairs  in  the  low  to  mid  latitude  ionosphere.  This  ionization  results  from  the  absorption 
by  neutrals  of  solar  extreme  ultraviolet  (EUV)  and  X-ray  radiation.  Production  due 
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Figure  43.  Major  and  minor  neutral  concentrations  during  solar  maximum  as  a  function 
of  altitude  as  computed  by  MSIS.  Heavier  ion  concentrations  decrease  with  altitude 
faster  than  the  light  ions 


to  photoionization  depends  on  the  incoming  solar  flux,  the  optical  depth,  the  neutral 
concentrations,  and  the  ionization  cross  section.  Each  of  these  factors  may  also  depend 
on  species,  altitude  (z),  wavelength  (A),  or  zenith  angle  (%);  for  a  particular  species 
s,  the  photoionization  rate  is  given  as 

r^si 

Ps(z,x)  =  Ns(z)  Co  (A)  exp  \—t(z,  x,  A)]  <x*(A)  dX  (20) 


84 


where 

/00(A)  is  the  solar  flux  at  the  top  of  the  atmosphere 
Ns(z)  is  the  number  density 
t(z,  x,  A)  is  the  optical  depth 

<7*  (A)  is  the  ionization  cross  section 

A si  is  the  threshold  wavelength  for  ionization 

Also  contributing  to  ion  production  are  chemical  reactions,  secondary  electron  pro¬ 
duction,  and  impact  ionization. 

The  primary  ion  loss  mechanisms  in  the  ionosphere  are  chemical  reactions.  The 
major  loss  reactions  include  ion-ion  recombination,  radiative  recombination,  dissocia¬ 
tive  recombination,  ion-atom  interchange,  and  associative  detachment.  An  example 
of  each  reaction  is  given  below. 


Ion-ion  Recombination: 

o2  +  o2  - 

— t  o2  +  o2 

Radiative  Recombination: 

0+  +  e~n  ■ 

— *  O  +  hv 

Dissociative  Recombination: 

N t  +  e~  - 

->•  N  +  N 

Ion-atom  Interchange: 

0+  +  H  - 

->  0  +  H+ 

Associative  Detachment: 

0~+0  — 

-t  02  +  e“ 

Production  and  loss  due  to  chemical  reactions  are  dependent  on  the  rate  of  the 
particular  reaction  and  on  the  concentrations  of  the  reactants.  For  a  reaction  such  as 

o+  +  n2^  n+  +  o 
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the  production  rate  of  iV^1"  (and  the  loss  rate  of  0+ )  is  given  as 


P(N+)  =  L{0+)  =  k  [0+]  [N2]  (21) 

where  k  is  the  reaction  rate  constant  for  the  specific  reaction  and  [0+]  and  [N2\  denote 
the  species’  concentration.  The  production  and  loss  rate  for  each  chemical  reaction  is 
the  terms  used  in  the  continuity  equation  (Equation  1).  The  reaction  rate  constant  k 
is  dependent  on  the  temperature  of  each  species  in  addition  to  the  activation  energy 
of  the  reaction,  which  can  be  described  as  the  minimum  energy  needed  to  form  the 
“new”  species  ( Schunk  and  Nagy,  2009).  The  reaction  rate  constant  is  determined  by 
the  equation 

kst  =  2d2^|^^^/ptexp(--^r\  (22) 

V  l‘st  V  hTstJ 

where  the  subscripts  denote  each  species  in  the  reaction,  d  describes  the  cross-section, 
kb  is  Boltzmann’s  constant,  fist  and  Tst  are  the  reduced  mass  and  temperature,  re¬ 
spectively,  and  Ea  is  the  activation  energy.  This  reaction  constant  varies  by  orders 
of  magnitude  for  the  different  reaction  types  listed  above  which  means  that  some 
processes  will  dominate  over  others.  For  example,  radiative  recombination  reaction 
rate  constants  are  on  the  order  of  10“ 12  while  dissociative  recombination  reaction  rate 
constants  are  on  the  order  of  10-7  ( Schunk  and  Nagy,  2009). 

In  the  D  region,  high  neutral  concentrations  lead  to  complex  chemical  reactions. 
One  of  these  reactions,  electron  attachment,  creates  the  only  significant  population  of 
negative  ions  in  the  ionosphere.  Deep  penetration  of  short  wavelength  X-ray  and  EUV 
radiation  also  contributes  to  the  complexity  of  the  D  region.  Due  to  this  complex 
chemical  nature  and  the  difficulty  in  modeling  the  shorter  wavelengths  of  the  solar 
spectrum,  the  D  region  is  not  included  in  the  Ionosphere  Forecast  Model  or  the 
Ionosphere-Plasmasphere  Model  and  will  not  be  discussed  further  here. 
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The  relatively  high  neutral  concentrations  in  the  E  region  lead  to  the  following 
chemical  reaction  production  mechanisms  for  NO+ 

N+  +  O  — >  NO+  +  N 
02  +  NO  — ►  N0+  +  02 
Ot  +  N2  — ►  NO+  +  NO 

Ion-atom  interchange  is  relatively  fast  compared  to  dissociative  recombination  in  this 
region  resulting  in  low  concentrations  of  0+  and  leaving  the  molecular  ions  (N2  , 
02,  and  NO+ )  as  the  major  ions.  In  the  Fi  region,  0+  becomes  an  important 
ion  due  to  the  higher  neutral  atomic  oxygen  concentrations  at  Fi  altitudes.  The 
reaction  rate  constants  determine  the  loss  mechanisms  in  these  regions.  Radiative 
recombination  is  slow  and  unimportant,  whereas  dissociative  recombination  occurs 
105  times  faster  ( Hargreaves ,  1992).  The  primary  dissociative  recombination  reactions 
in  the  E  and  Fi  regions  are 


02  +  e  — 

O  +  O 

Nt  +  e~  - 

-+  N  +  N 

NO+  +  e~  - 

-+  N  +  O 

The  F2  region  has  the  same  production  and  loss  processes  as  the  E  and  Fi  regions, 
however,  at  and  above  F2  region  altitudes,  plasma  transport  is  significant  and  deter¬ 
mines  the  density  structure.  While  the  F2  region  is  controlled  by  all  three  processes, 
the  topside  ionosphere  is  controlled  by  transport  alone. 


87 


2.3  Transport 


Transport  in  the  ionosphere  is  governed  by  the  momentum  equation 


As" 


Dsus 

~dT 


—  —  Vps  —  V  •  ts  +  psG  +  nsqs 


E 


us  x  B 


^  '  Ps^st  Us)  (23) 


where  the  subscripts  denote  the  species  s  and  the  target  t,  ps  is  the  mass  density, 
Ds/Dt  is  the  convective  derivative,  ps  is  the  pressure,  kb  is  Boltzmann’s  constant,  ts 
is  the  stress,  and  vst  is  the  collision  frequency.  In  obtaining  Equation  23,  thermal 
diffusion,  Coriolis  force,  and  centripetal  force  were  all  assumed  to  be  negligible.  After 
also  making  the  diffusion  approximation  (steady  state  and  subsonic)  Equation  23 
reduces  to 


Vps  +  V  •  rs  -  psG  -  nsqs 


E  +  us  x  B 


y.  ps^st  ( ut 


U  a 


(24) 


In  order  to  further  simplify  the  solution  to  the  momentum  equation,  two  separate  cases 
are  considered.  First,  diffusion  along  the  magnetic  field  will  be  examined  followed 
by  diffusion  perpendicular  to  the  magnetic  field.  The  overall  plasma  motion  is  a 
combination  of  both  solutions. 

For  diffusion  along  the  magnetic  field,  the  plasma  is  confined  to  orbit  the  magnetic 
field  lines.  In  this  situation  the  electrons  and  the  major  ion  species  move  together  and 
charge  neutrality  (ne  =  nt)  and  zero  current  (neue  =  npu^)  conditions  prevail.  This 
ambipolar  diffusion  behavior  is  a  result  of  the  polarization  electric  field  that  develops 
due  to  the  slight  charge  separation  of  the  electrons  and  ions  due  to  their  different 
masses.  With  this  ambipolar  diffusion  approximation  and  the  previously  mentioned 
assumptions,  the  ion  and  electron  momentum  equations  along  the  magnetic  field  are 


given  as 


V||  Pi  H-  Pi& ||  Pi^ie  (^e  ^i)  ||  “1“  Pi^in  {^n  ^i)  ||  (25) 

V|| Pe  A  *  ^~e ) 1 1  Pe^\\  H-  ^e^-^||  Pe^ei  ^e)||  “I”  Pe^en  ip^n  ^e)||  (25) 


where  E\\  is  the  polarization  electric  field.  Adding  these  equations  together  and 
solving  for  the  ion  velocity  gives  the  ambipolar  diffusion  equation 


^i\\  Da 


— Vnrq  + 

Hi 


_ VmT 

rji  V| \1P 

J.  n 


rnjG_ ||  (V  '  T») || 
2kf,Tp  2riikbTp 


m 

2kbTp 


^in^nl  I 


(27) 


where  terms  containing  the  electron  mass  (me)  were  neglected  and  the  ideal  gas  law 
was  used,  and  where  the  ambipolar  diffusion  coefficient  (Da)  and  plasma  temperature 
(Tp)  are  given  by 


Da 

TP 


2kbDp 

Te  +  T\ 

2 


Equation  27  shows  that  in  steady  state,  the  ions  will  move  along  the  magnetic  held 
with  the  neutral  wind  subject  to  temperature  gradients,  density  gradients,  gravity, 
and  stress  forces. 

For  diffusion  across  the  magnetic  held,  the  stress  term  in  Equation  24  is  also 
assumed  to  be  negligible  giving 


Vps  -  psG  -  nsqs 


E_ i_  us  x  B 


Ps^sn  (^n  Us) 


(28) 


where  E±  is  an  applied  electric  held,  the  subscript  s  denotes  either  an  electron  or 
ion,  and  the  electron-ion  collisions  are  neglected  because  the  momentum  transfer 
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is  small  ( Schunk  and  Nagy ,  2009).  Transforming  Equation  28  to  a  reference  frame 
moving  with  the  neutral  wind  (us  — >•  us  +  un)  simplifies  the  solution  and  introduces 
an  effective  electric  field  given  by  E'±  =  E±  +  un  x  B.  Equation  28  then  becomes 


Vps  -  psG  -  nsqs 


E  i  +  us±x  B 


Ps^sn^sA. 


(29) 


Solving  for  usJ_  gives 


u 


si  I  ^s  I 

msvsn  ' 


x  B  = 


-1 


~V±Ps  +  —G±  +  Z^E± 


TAIsL'sn 


(30) 


Equation  30  is  simplified  as 

K±  -  K  (K±  x  b)  =  — —V±ps  +  —G±  +  psE'±  (31) 

'  Ps  V sn 

using  the  following  relations 


Ps 

nskbTs 

ps 

Qs 

_  qsB 

B 

cs 

ms 

b  = 

B 

Ds 

kbTs 

K  = 

^ cs 

Wls^sn 

V sn 

Because  usJ_  is  perpendicular  to  the  magnetic  field,  it  can  have  two  orientations. 
Consider  a  coordinate  system  at  the  equator  where  the  magnetic  field  B  is  pointed 
north  in  the  b  direction.  Then  the  components  of  us±  can  be  in  the  downward  (i) 
direction  or  the  horizontal  eastward  (— j)  direction.  Separating  us±  in  Equation  31 
into  these  components  and  simplifying  gives 

~D  .  1  ^ 

i  — - V ±Dpsi  H  G . 

Ps  vsn 


K±P  -  K±J  +  kK±J  +  kK±h 


Di  +  psE±Di  (32) 


90 


where  the  variations  of  the  density  gradient,  gravity,  and  electric  field  are  small  in  the 
horizontal  direction  and  can  be  neglected.  Equating  all  the  terms  in  the  i  direction 
and  all  the  terms  in  the  j  direction,  we  can  write 


- V ±ps  H - G±_  +  psE j 

Ps  E sn 


(33) 

(34) 


Substituting  us±H  from  Equation  34  into  Equation  33  and  simplifying  gives 

US±D  =  — y  +  — G±  +  nsE_ L  (35) 

Equation  35  describes  the  vertical  drift  that  is  perpendicular  to  the  magnetic  field 
but  parallel  to  the  force  causing  the  drift.  For  G±  =  E'±  =  0  in  an  unmagnetized  case 
(K  =  0),  Equation  35  reduces  to  Fick’s  Law,  which  states  that  the  particle  flux  is 
proportional  to  the  density  gradient.  Substitution  of  Fick’s  Law  into  the  continuity 
equation  results  in  the  classical  diffusion  equation  ( Schunk  and  Nagy ,  2009). 

The  horizontal  drift  is  obtained  by  substituting  Equation  35  for  us±D  back  into 
Equation  34  and  simplifying,  which  results  in 


Us±„  =  l  +  yK 2  [«P  +  +  UE] 


where 


-1  (V-LPS  x  B\ 

up  = - - - - - —  is  the  gradient  drift 

nsqs  B2 

ms  x  B) 

ug  =  — — — — —  is  the  gravitational  drift 

Qs  B2 

(El  x  B) 

up  =  — — — — —  is  the  electromagnetic  drift 
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These  three  drifts  are  perpendicular  to  both  the  magnetic  field  and  the  force  (pressure 
gradient,  gravity,  or  the  electric  field)  causing  the  drift.  Electrons  and  ions  drift 
across  the  magnetic  field  in  opposite  directions  in  the  presence  of  pressure  gradients 
and  gravity,  but  they  drift  together  in  the  presence  of  a  perpendicular  electric  field. 
Typically  the  drifts  due  to  pressure  gradients  and  gravity  are  small  and  the  dominant 
drift  is  the  electromagnetic  drift. 

Since  usJ_  =  u's±D  +  us±  ,  the  total  plasma  drift  across  the  magnetic  field  is 


1 

1  +  A'2 


- Vj _ps  H - G±  +  ^SE_ |_ 

Ps  V  sn 


1 


1  + 


1 

K2 


[Up  +  Uq  +  Up\ 


(37) 


When  collisions  are  negligible  K  —>  oo  and  the  first  term  of  Equation  37  drops  out 
leaving  only  the  motion  perpendicular  to  both  the  magnetic  field  and  the  force.  This 
approximation  is  valid  at  high  altitudes.  In  locations  where  K  — >•  0,  the  second  term 
drops  out  leaving  only  the  motion  perpendicular  to  the  magnetic  field  but  parallel  to 
the  force.  This  situation  occurs  at  low  altitudes  where  the  magnetic  field  strength  is 
small  compared  to  the  collision  frequency.  At  intermediate  locations  where  neither 
approximation  is  valid,  the  plasma  motion  will  be  a  combination  of  these  motions  and 
will  move  perpendicular  to  the  magnetic  field  and  at  an  angle  to  the  force  causing 
the  motion  ( Kelley ,  1989). 
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